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ABSTRACT 

Reliably  modeling  noise  attenuation  through  interaction  with  vibrating  bound¬ 
ary  structures  is  fundamental  to  the  formulation  of  effective  active  noise  control  sys¬ 
tems.  In  this  paper  we  investigate,  through  numerical  approximation,  uniform  ex¬ 
ponential  stability  of  two  systems  which  model  the  acoustic/structure  interaction  of 
an  air-filled,  rectangular  cavity.  The  first  model  assumes  dissipative  boundary  condi¬ 
tions  along  one  side  of  the  boundary,  while  the  second  assumes  dissipative  boundary 
conditions  along  all  four  sides  of  the  cavity.  We  obtain  weak  variational  formulations 
for  these  models,  express  each  as  finite  dimensional  systems,  and  use  the  Galerkin 
technique  to  transform  the  distributed  parameter  systems  into  systems  of  ordinary 
differential  equations.  We  analyze  the  stability  of  the  finite  dimensional  systems  in 
order  to  gain  insight  into  the  stability  of  the  original  infinite  dimensional  systems. 
Essentially,  our  analysis  consists  of  solving  a  generalized  eigenvalue  problem  and  ob¬ 
serving  where  the  eigenvalues  lie  within  the  complex  plane.  This  stability  analysis 
leads  us  to  conclude  that  one  model  is  better  suited  for  use  in  the  formulation  of  the 
noise  control  problem 


V 


TABLE  OF  CONTENTS 


I.  INTRODUCTION .  1 

II.  MATHEMATICAL  MODELS  . 5 

A.  MATHEMATICAL  MODEL  I .  8 

B.  MATHEMATICAL  MODEL  II  .  16 

III.  GALERKIN  APPROXIMATION  METHOD .  21 

IV.  FINITE  DIMENSIONAL  APPROXIMATIONS .  29 

V.  SPECIFIC  APPROXIMATIONS  AND  RESULTS .  35 

VI.  CONCLUSIONS .  53 

APPENDIX.  MATLAB  FUNCTION  AND  SCRIPT  FILES .  55 

REFERENCES .  75 

INITIAL  DISTRIBUTION  LIST  .  77 


LIST  OF  TABLES 


I.  Model  I:  Margin  between  the  open  loop  eigenvalues  (A)  and 

the  imaginary  axis .  46 

II.  Model  II:  Margin  between  the  open  loop  eigenvalues  (A)  and 

the  imaginary  axis .  46 


IX 


X 


LIST  OF  FIGURES 


1.  2-D  Acoustic  Chamber  for  the  General  Mathematical  Model  .  .  5 

2.  2-D  Acoustic  Chamber  for  Model  I .  8 

3.  2-D  Acoustic  Chamber  for  Model  II .  16 

4.  Cubic  Polynomial  Splines .  37 

5.  Transformed  Legendre  Polynomials .  38 

6.  Mod  I:  Eigenvalues  for  n  =  rtix  =  my  =  6,7, 8, 9, 10, 11, 12, 13.  .  .  47 

7.  Mod  I:  Eigenvalues  for  n  =  =  my  =  lA,lh, 16 .  48 

8.  Mod  I:  Eigenvalues  for  n  =  m^  =  my  =  17, 18 .  49 

9.  Mod  II:  Eigenvalues  for  n  =  m^  =  m.y  =  5, 6, 7, 8, 9, 10, 11, 12.  .  .  50 

10.  Mod  II:  Eigenvalues  for  n  =  m^  =  my  =  13, 14, 15, 16, 17, 18, 19, 20.  51 


XI 


I. 


INTRODUCTION 


Imagine  that  you  are  ten  years  old  and  your  teenage  brother  asks  you  to 
participate  in  a  scientific  experiment.  You  agree.  Inside  a  55  gallon  drum  you  go. 
Your  brother,  clever  lad  that  he  is,  manages  to  suspend  the  drum  (and  you)  a  few  feet 
in  the  air.  So  there  you  are,  dangling  by  a  rope  three  feet  off  the  ground-sealed  in  a 
metal  drum.  With  your  baseball  bat  in  hand,  he  strikes  the  drum  dead  center.  After 
the  drum  stops  reverberating,  he  lowers  it  to  the  ground.  As  he  pops  the  lid  off,  he 
asks  you  to  describe  how  the  acoustic  noise  field  in  the  drum  behaved  after  his  forceful 
swing.  As  a  composed  budding  young  scientist,  you  disregard  the  blood  dripping  from 
your  left  earlobe  and  answer,  “Initially,  a  large  noise  field  was  created  by  the  strike  of 
the  bat.  However,  after  the  impulse  force  was  experienced,  the  intensity  of  the  noise 
field  steadily  decreased  until  it  eventually  was  imperceptible.”  Your  brother  thanks 
you  dearly  and  then  proposes  a  heat  conduction  experiment... 

In  this  paper,  we  too  examine  noise  transmission  attenuation  through  a  vi¬ 
brating  boundary  structure.  Our  approach,  however,  is  to  do  so  through  numerical 
approximation.  More  specifically,  we  examine  the  exponential  stability  of  infinite  di¬ 
mensional  second  order  systems  of  coupled  partial  differential  equations  by  numerical 
approximation.  This  is  a  topic  attracting  considerable  interest  in  the  fields  of  engi¬ 
neering  and  applied  mathematics  because  the  applications  are  both  numerous  and 
diverse — reducing  noise  levels  in  automobiles,  aircraft,  and  space  launch  vehicles  to 
name  a  few.  The  degree  to  which  a  control  system  formulated  to  effect  noise  reduc¬ 
tion  in  a  fluid-filled  cavity  succeeds  depends  to  great  extent  on  how  accurately  the 
underlying  mathematical  model  agrees  with  the  observed  physical  phenomena.  As 
our  young  scientist  reported  above,  the  acoustic  field  created  by  an  external  force 
acting  on  the  cavity  boundary  steadily  diminished  to  zero  as  time  passed  (in  the 
absence  of  any  sustaining  force).  This  is  equivalent  to  requiring  that  all  solutions  to 
the  system  of  equations  selected  to  model  the  behavior  of  the  acoustic  field,  as  well 
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as  the  vibrating  boundary,  converge  exponentially  to  zero  by  a  uniform  rate  of  decay. 
Establishing  this  result  for  the  infinite  dimensional  system  is  nontrivial  and,  at  times, 
a  very  difficult  task. 

To  help  clarify  these  ideas,  consider  the  following  abstract  formulation.  Ac¬ 
cording  to  [Ref.  1],  many  examples  related  to  acoustics  or  fluid/structure  interactions 
can  be  modeled  abstractly  as  a  first  order  system  of  equations  as  follows: 

yt{t)  =  Ay{t),  t>0,  j/(t)e^,  (I.l) 

where  H  is  an  appropriately  defined  Hilbert  space.  The  companion  linear  control 
system  is  typically  written 

y,{t)  =  Ay{t)  +  Bh{t),  h{t)e'R\  (1.2) 

where  /i  is  a  control  input  and  R  is  a  linear  operator  from  'RA  into  H.  This  system 
possesses  uniform  exponential  stability  if  there  exist  M  >  0  and  >  0  such  that  for 
alH  >  0  and  for  all  (j/(0),  ?/t(0))  €  % 

<  -^e"^*||(y(0),t/t(0))||w  [Ref. 2,  3] , 

where  \\y(t),yt{t)\\-u  denotes  the  energy  of  the  system  at  time  t  and  ||(y(0), yi(0)||7^ 
denotes  the  energy  of  the  system  at  t  =  0. 

The  authors  of  [Ref.  1]  state  that  “the  most  common  approach  for  the  ap¬ 
proximation  of  a  control  problem  involving  1.2  is  to  formulate  a  sequence  of  finite 
dimensional  control  systems  of  the  form 

yfit)  =  A^y^it)  +  B^h{t) ,  f  >  0 ,  y^ {t)  G  ,  (1.3) 

where  the  dimension  of  the  finite  solution  space  increases  toward  infinity  cis  N 
tends  to  infinity.  In  general,  equation  1.3  is  derived  from  1.2  using  space  discretization 
techniques  such  ais  finite  difference,  finite  elements  or  spectral  methods  developed 
for  the  approximation  of  the  solutions  of  I.l.  A  control  strategy  is  then  designed 
for  the  finite  dimensional  control  problem  involving  1.3.  This  control  is  used  as  an 
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approximation  to  the  desired  control  function  for  the  infinite  dimensional  control 
problem  1.2.  One  of  the  most  practical  conditions  to  assure  the  well-posedness  of 
the  finite  dimensional  control  problem,  as  well  as  the  convergence  of  the  approximate 
controls  to  the  desired  control  for  the  infinite-dimensional  system,  is  that  the  solutions 
of  1.3  for  =  0  preserve  the  exponential  decay  of  the  solutions  of  I.l.”  Hence 
determining  the  stability  of  the  finite  dimensional  system 

y^{t)=z  t>0,  (1.4) 

is  of  fundamental  importance. 

Typically  this  stability  analysis  is  accomplished  by  examining  where  the  eigen¬ 
values  of  he  within  the  complex  plane,  since  theory  tells  us  that  system  1.4  is 
globally  stable  if  and  only  if  all  eigenvalues  of  A^  have  negative  real  parts  [Ref.  4]. 
The  process  is  nontrivial  because  numerical  results  have  indicated  that  many  popular 
approximation  schemes  fail  to  maintain  a  uniform  decay  rate  as  the  dimension  of  the 
approximating  system  1.4  increases  [Ref.  1],  even  in  cases  where  the  original  system 
was  proven  to  be  exponentially  stable. 

In  this  paper,  we  examine  two  different  infinite  dimensional  models  by  numer¬ 
ical  approximation  in  order  to  gain  insight  into  their  adequacy  for  control  system 
formulations.  Model  I  is  believed  to  be  stable,  but  not  uniformly  exponentially  sta¬ 
ble  [Ref.  5].  Model  II  is  believed  to  be  exponentially  stable  [Ref.  6]. 

In  the  following  chapters,  topics  introduced  above  are  addressed  in  greater 
detail.  In  Chapter  II,  we  develop  a  general  and  two  specific  models  describing  an 
acoustic  field  inside  a  two-dimensional  fluid-filled  cavity  surrounded  by  a  perturbable 
boundary.  In  Chapter  III,  we  illustrate  the  Galerkin  technique  chosen  for  our  numer¬ 
ical  approximations,  and  in  Chapter  IV,  we  apply  this  technique  to  the  models  under 
consideration.  In  Chapter  V,  we  describe  how  specific  approximations  were  obtained 
and  present  our  results.  We  conclude  with  summary  comments  and  propose  future 
areas  of  study  in  Chapter  VI. 
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II. 


MATHEMATICAL  MODELS 


Force 


r=i;ui^ui; 


Figure  1.  2-D  Acoustic  Chamber  for  the  General  Mathematical  Model 


In  this  chapter,  we  develop  a  general  and  two  specific  models  for  a  system  con¬ 
sisting  of  the  wave  equation  coupled  with  the  beam  equation  on  part  of  the  boundary. 
Consider  the  two-dimensional  rectangular  air  filled  cavity  surrounded  by  an  impene¬ 
trable  boundary  shown  in  Figure  1.  A  noise  source  exterior  to  the  cavity  produces  a 
perturbing  force  /  which  induces  vibrations  in  the  cavity  boundary  causing  fluctua¬ 
tions  (i.e.,  undesirable  noise)  in  the  acoustic  pressure  field  within  the  cavity.  Equa¬ 
tions  of  motion  describing  boundary  vibrations  and  acoustic  pressure  fluctuations 
within  the  cavity,  together  with  appropriate  initial  value  and  boundary  conditions 
form  a  system  of  coupled,  second  order  partial  differential  equations  in  time. 

We  begin  formulation  of  the  general  model  by  assuming  that  the  interior  acous¬ 
tic  pressure  field  satisfies  the  standard  wave  equation  <f>tt  =  c^A(f>  where  is  the 
velocity  potential  throughout  the  cavity  and  c  is  the  uniform  speed  of  sound  in  the 


fluid.  The  velocity  potential  ^  is  a  complex-valued  function  satisfying 

v{t,x,y)  = 

where  v  denotes  the  fluid  velocity.  Initial  value  conditions  x,  y)  and  x,  y),  as 
well  as  boundary  conditions  along  dCl — either  Dirichlet,  Neumann  or  dissipative — are 
specified. 

A  Newtonian  analysis  of  forces  and  bending  moments  leads  to  equations  de¬ 
scribing  the  motion  of  the  elastic  walls  bounding  the  cavity.  For  simplicity  we  assume 
that  the  boundary  walls  are  impenetrable  and  that  only  one  side  of  the  rectangular 
boundary  is  perturbable  (The  term  beam  refers  to  the  perturbable  boundary.).  We 
assume  an  Euler-Bernoulli  beam  where  (i)  w{t,  x)  denotes  the  transverse  displace¬ 
ment  of  the  beam  of  length  a,  (ii)  pt,  and  p/  denote  the  uniform  mass  densities  of  the 
beam  and  fluid,  respectively,  (iii)  M{t,  x)  denotes  the  total  internal  moment  of  the 
beam,  and  (iv)  f{t,x)  denotes  the  external  forcing  term.  The  beam  equation  takes 
the  general  form 

pf,Wttit,x)  +  Mxxit.x)  =  -pf(j)t{t,x,w{t,x))  +  f{t,x)  ,  (II.l) 

where  —pf(l)t{t,x,w{t,x))  is  the  acoustic  pressure  (This  is  the  first  coupling  term  we 
see  in  the  acoustic/structure  system.). 

Initial  value  conditions  are  specified,  and  boundary  conditions  for  the  beam 
indicate  whether  the  ends  are  free,  partially  restrained,  or  clamped.  Additionally, 
assumptions  specifying  the  types  of  internal  moments  of  the  beam  are  necessary. 
Internal  moments  typically  consist  of  bending  and  damping  moments, 

M (t,  x)  =  Mbending  T  ^damping  5 

where 

Mbending  =  strain  Component  =  E{x)I{x)wxx{t,x) . 

The  stiffness  of  the  beam  is  given  by  E{x)l{x),  where  E{x)  is  the  Young’s  modulus 
and  /(x)  is  the  cross-sectional  area  of  the  beam.  Mdamping  is  taken  to  be  either 
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Kelvin- Voigt,  spatial  hysteresis,  or  time  hysteresis  damping.  While  spatial  hysteresis 
damping  is  an  appropriate  choice  for  a  beam  constructed  of  composite  materials  and 
time  hysteresis  damping  is  an  appropriate  choice  for  a  beam  constructed  of  “material 
with  memory” ,  Kelvin- Voigt  damping  assumes  a  memoryless  beam  of  uniform  (linear) 
mass  density  where  the  damping  stress  is  proportional  to  the  strain  rate.  That  is, 

Kelvin- Voigt  damping:  Mdamping  =  strain  rate  component  =  CD{x)I{x)wxxt , 

where  C£)/(x)  is  the  product  of  the  Kelvin- Voigt  damping  coefficient  Cd{x)  and  beam 
cross-sectional  area  I{x)  [Ref.  7].  Here  we  consider  a  beam  of  uniform  cross  sec¬ 
tional  area  and  density.  Thus  we  assume  Kelvin- Voigt  damping  and  take  the  coeffi¬ 
cient  functions  cd{x)I{x)  and  E{x)I{x)  to  be  constant  for  all  x  (i.e.,  cb{x)I{x)  — ^ 
C£)I  and  E{x)I{x)  — v  El). 

Based  upon  the  above  discussion,  the  generalized  mathematical  model  for  the 
acoustic/structure  system  is: 

^tt  =  for  {x,y)  G  H  t  >  0, 

(II.2) 

PbWu{t,x)  -I-  Mxx{t,x)  = -pf4>t{t,x,w{t,x))-\- f{t,x) ,  d  <  X  <  a ,  t>0, 

with  appropriate  initial  value  and  boundary  conditions  specified  for  the  wave  and 
beam  equations. 

Next  we  present  two  specific  models  and  formally  show  that  each  model  can 
be  expressed  in  both  weak  and  strong  formulations — which,  as  we  shall  see,  are  equiv¬ 
alent  given  appropriate  choices  of  inner  product  spaces.  This  preliminary  work  will 
formally  justify  expressing  Models  I  and  II  as  first  order  systems  in  time  thus  facil¬ 
itating  our  numerical  examinations  of  solution  stability.  In  Chapter  III  we  examine 
in  some  detail  the  popular  variational  scheme,  the  Galerkin  method,  used  to  obtain 
approximate  solutions  to  coupled  systems  such  as  those  addressed  in  this  paper.  We 
use  the  Galerkin  scheme  to  transform  the  infinite  dimensional  system  II.2  into  a  finite 
dimensional  system  to  facilitate  numerical  analysis. 
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A.  MATHEMATICAL  MODEL  I 


Force 

Figure  2.  2-D  Acoustic  Chamber  for  Model  I 

The  theoretical  model  considered  in  this  section  is  shown  in  Figure  2.  Here  the 
wave  equation  is  coupled  with  the  beam  equation  on  one  side  of  a  2-D  rectangular, 
air  filled  cavity.  We  assume  Neumann  boundary  conditions  along  F.  That  is, 

V(j>  •  n  =  0,  (x,y)  e  r,  t  >  0, 

where  n  represents  the  the  outer  normal  with  respect  to  F.  The  boundary  is,  in 
effect,  a  sound-insulated  rigid  wall  which  prevents  any  acoustic  energy  from  escaping, 
or  alternatively,  the  boundary  is  a  perfect  reflector  of  acoustic  waves. 

Further,  we  assume  that  the  perturbable  boundary  (i.e.,  beam)  can  be  charac¬ 
terized  as  an  impenetrable  fixed-end  Euler-Bernoulli  beam  with  Kelvin- Voigt  damping 
and  that  both  ends  of  the  beam  are  clamped.  Given  these  assumptions,  the  equations 
of  motion  describing  the  vibrations  of  the  perturbable  boundary  are: 

phWtt{t,x) =  -pj(f)t{t,x,w{t,x))  +  f{t,x) ,  0<x<a,  t>0  (II.3) 
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u;(i,0)  =  0)  —  w{t,a)  =  Wx{t,a)  =  0,  t  >  0, 


(11.4) 


where 


'w{t,x) 


Pb 


PS 

M{t,  x) 


the  transverse  displacement  of  the  beam 
the  fluid  velocity  potential 
the  linear  mass  density  of  the  beam 
the  uniform  mass  density  of  the  fluid 
El XI) XX  “h  X)j^Ixi)xxt 

the  force  due  to  an  exterior  noise  field. 


For  the  stability  analysis  of  this  model,  we  consider  the  open  loop  problem  absent  any 
exterior  noise  field  (i.e.,  f{t,x)  =  0  for  t  >  0)  and  do  not  concern  ourselves  with  any 
noise  control  aspects.  We  assume  (i)  that  the  beam  is  impenetrable  to  the  adjoining 
fluid  and  obtain  the  second  coupling  equation  (i.e.,  the  continuity  of  velocity) 

'U)t(x,t)  =  V(j>{t,x,xv{x,t))  -  n,  0<a:<a,t>0,  (II.5) 

and  (ii)  that  displacements  from  the  beam’s  position  of  rest  are  small,  which  is  in¬ 
herent  in  the  Euler-Bernoulli  formulation.  Because  of  (ii),  we  take  the  transverse 
displacement  of  the  beam  as  w(t,x)  =  w(t,x)  +  S  where  w  =  0.  Under  these 
assumptions,  equation  II.3  and  II.5  become 

pi,Wtt(t,x)  +  Mxx(t,x)  =  -Pf(^t(t,x,w(t,x)  +  S)  (II.6) 

wt(t,x)  =  V(j){t,x,w{t,x)  -h  d)  •  n  .  (11.7) 

By  using  two  term  Taylor  Series  expansions  of  (fxt  and  with  respect  to  y  about  the 
point  X,  equations  II.6  and  II.  7  become 

PbV)u{t,x)  -1-  Mxx{t,x)  =  -Pf[(f)t{t,x,0)  +  (f)tyit,X,0)w], 
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Wt{t,x)  =  V<^(t,a:,0)  •  n  +  {'V4>y{t,x,0)w)  ■  n. 

We  drop  the  higher  order  terms  —pf(f)ty[x^Q,t)w  and  (V<^y(a;,  0,  •  n  in  these 

two  equations  because  of  the  assumption  of  small  beam  displacement  and  obtain  first 
order  approximations  for  —pf<j>t  and  V<f> ,  respectively.  Upon  approximating  the  space 
domain  Q(t)  by  fi  =  [0,  a]  x  [0,^],  the  open  loop  model  described  above  is  given  by 

(j)tt  =  {x,y)  e 

V(j)  ■  n  =  0  ,  (a;,  ?/)  €  r  ,  t  >  0, 

(l>y{t,x,0)  =  -wt{t,x),  0<x<a,t>0, 

PiWtt  +  dl{EIWxx  +  CDiWxxt)  =  -pf<l>t{t,x,0) ,  0  <  X  <  a,  t  >  0,  /  (H-S) 

■U7(t, 0)  =  0)  =  w{t,a)  =  Wx{t,a)  =  0,  t  >  0, 

4>{0,x,y)  =  <f>o{x,y)  ,  w{0,x)  =  wo{x), 

<f)t{0,x,y)  =  (f>i{x,y)  ,  wt(0,a;)  = 

Note:  Throughout  this  paper,  da  denotes  partial  dilferentiation  with  respect  to  the 
variable  a  (e.g.,  = 

System  II.8  is  a  formal  representation  of  the  dynamics  of  a  coupled  acous- 
tic/beam  structure.  Computational  techniques  (e.g.,  variational  methods)  used  to 
obtain  approximate  solutions  to  this  system  are  based  on  rigorous  convergence  argu¬ 
ments,  typically  done  in  the  context  of  variational  formulations  of  II. 8.  To  accomplish 
this,  the  state  is  taken  to  be  z{t)  —  (4>,  w)  in  the  Hilbert  space  H  =  L  (H)  x  L^{To) 
with  energy  product 

where  L  (H)  is  the  quotient  space  of  over  the  constant  functions.  Also,  we  define 
the  Hilbert  space  V  =  x  HQ(ro)  where  is  the  quotient  space  of 

over  the  constant  functions  and  Ho(ro)  =  {V’  €  H^(ro)  :  V’(x)  =  i^x{^)  =  0  at 
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X  =  0,  a  }.  The  energy  product  of  V  is  taken  as 


\Wy 


?  I  I  )  /  p ^  ^ ^ ^xx'pxx  ^'y  • 

\ri  j  /  ^  J^o 

Next  we  define  the  weak  variational  (i.e.,  sesqnilinear)  forms: 

=  ^P/ W- for  <f),^  eTT^iO.) 

/32{w,rj)=  /  Wa^xVxxd'y  for  w,r]  e  HoiTo) 

''To 

pi{w,T})=  /  pbwrjd'y  for  w,t)  e  L^iTo) 

JTq 

P2{<l>,0=  f  do;  for  ^,^6  1^(0) 

pi{w,ri)=  /  Elw^^rj^^d'y  for  10,77  e /^^(ro) 

•/To 

A‘2(^, 0  =  j^Ps^(t>  ■'^idio  for  (/), ^  e  ^\f]) 

«i(u;,?7)=  /  CDlWxx'Qxxd'^i  for  w,r}  E  Hl{Vo) 

^)  =  /  P/  0)?7C?7  for  (j)  6  F^(n)  and  77  €  /fo(ro) 

''Fq 

'r2(«^,0=  /  2:5  0)to  ^7  for  to  G /fo(ro)  and  ^ 

-'Fo 


and  express  system  11.8  in  weak  form  as, 


(II.9) 


(11.10) 

(11.11) 


Our  next  task  is  to  write  equations  11.10  and  II.  11  as  a  single  second  order 
differential  equation.  Let  =  (<7;>,u;)  and  ^  =  (^,77),  such  that  G  V,  and  define 
sesquilinear  forms: 

<7i($,^)  =  P2+M1  =  J^PfV<^-V^du)  +  EIw^x'Hxxd'i , 

o’2(^,^)  =  Ki  +  ri-r2  =  f  {cdIwxxVxx  +  0)77  -  ^(t,  x,  0)u;}  ^7  , 

«^Fo 

where  ai,(T2  G  V  x  V  — y  C  (space  of  complex  numbers).  The  formulations  ai  and 
£72  satisfy  coercivity  and  continuity  (i.e.,  boundedness)  conditions 

5Rcri($,$)  >  cill^ll^. 
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Ik.($.-l>)||  <  C2||#i|.,||<P||v. 

>  C3{TO,r,  TO„)iJ(r„)  =  CsIIuiIIhJP,)  , 

Il<r2($,3')i|  <  C4||$||v||'I'|k, 

where  3?  denotes  the  “real  part  of” .  The  second  order  open  loop  problem  is  given  by 


{ztt{t),  ^)y.,y  +  cr2{zt{t),  ^)  =  0  ,  (11.12) 

where  V*  is  the  dual  of  space  V  and  (•,  •)v,v  denotes  the  usual  duality  pairing. 

Since  <Ti  and  (72  satisfy  the  continuity  and  coercivity  conditions  shown  above, 
the  Lax-Milgram  Theorem  guarantees  the  existence  of  uniquely  determined  bounded 
linear  operators  Ai,A2  such  that  the  weak  and  strong  formulations  of  the  coupled 
system  are  equivalent  [Ref.  8].  That  is 


(Ai$,  =  (7i($, '1?)  and  (^4.2$,  =  o’2($,  ^) . 


Thus  system  11.12  gives  rise  to  the  equivalent  system  defined  in  terms  of  functionals 
Ai  and  A2 

+  A2Zt{t)  +  Axzit)  —  0 .  (11.13) 


To  facilitate  our  numerical  work,  we  must  express  system  II. 8  as  a  first  order 


1  is  to  write  Model  I 

as 

4>t 

0 

0 

I 

0 

Wt 

0 

0 

0 

I 

4>tt 

c^A 

0 

0 

0 

.  . 

0  - 

-—dt 

Pb  ^ 

-n 

Pb  ^ 

<t> 

W 

.  . 

(11.14) 


u,  A 

The  symbol  11  appearing  in  matrix  A  above  represents  whenever  the 

product  Au  is  calculated. 

We  now  develop  the  necessary  notation  and  sesquilinear  forms  in  order  to  write 
system  II. 8  as  a  first  order  system.  The  following  approach  replicates  that  found  in 
[Ref.  9]. 
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We  begin  by  defining  the  product  spaces  V  =  V  xV  and  H  =  V  x  H  whose 
norms  are  given  by 

11(^5  ^)|lv  =  ll^lly  +  ll^llv  and  j|($,  ^)||^  =  ||$||^  +  ||^|[^  ^  respectively. 

For  X  =  (^, '!')  and  0  =  (T,  A),  the  sesquilinear  form  cr  :  V  x  V  — y  C  is  then  defined 
by 

a(0,x)  =  a((T,A),($,^))  =  -(A,$)v  +  o-i(T,^)  +  o-2(A,^). 

Since  the  duality  product  (•,  ■)v,v  is  the  unique  extension  by  continuity  of  the  scalar 
product  (•,  ■)ff  from  H  x  V  to  V*  x  F,  it  follows  that  for  appropriate  restrictions  on 
0  we  can  write 

^(0,x)  =  c^((T,A),($,^))  =  -(A,#)y  +  (AiT,^)v.,v  +  (A2A,^)v.,v 

=  — (A,  $)y  +  (AiT  +  A2A, 

=  ((—A,  AiT  +  A2A),  ($, '!'))•« 

=  {-Ae,xU. 

The  operator  A  :  %  — >  'H  is  given  by 

A  = 

where  the  domain  of  A  =  {0  =  (T,  A)  G  :  A  e  F,  AiT  +  A2A  G  H},  Ai  and  A2 
are  the  operators  defined  by  <ti  and  (T2,  respectively,  and  the  above  calculations  hold 
for  0  G  domain  A. 

By  letting  Z{t)  =  (z{t),zt{t))  and  taking  x  €  V,  the  first  order  system  is 
written  in  weak  form  as 


0  I 
-Ai  — A2 


(11.15) 


{^tit)^x)v-,v  =  -a{Z{t),x), 


which  is  formally  equivalent  to  the  strong  formulation 


Zt{t)  =  AZ{t) 


(11.16) 
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in  1-L,  where  A  is  given  in  11.15.  Equation  11.16  concisely  represents  the  matrix 
representation  of  Model  I  given  by  equation  11.14  where 


A  = 


0  0/0 
0  0  0  / 

c^A  0  0  0 

0  -n 

Ph  ^  Pb  ^ 


and 


-c^A  0 

0 

0 

Ai  — 

,  A2  — 

i 

O 

n 

cnlai 

Pb  ®  J 

The  linear  operator  A  is  dissipative  in  the  sense  {Axi  x)'H  <  0  for  x  G  domain  A.  To 
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see  this  let  x  =  (^,  w,  (f)t,  Wt),  then 


{•^X,  x)n  = 


4>t 

Wt 

c^A<f> 

f  9>  -  n^,  - 


’  ^  ' 

w 

-  _ 

n^vxH 


=  </>/V0(,  V</>)x,2(n)  +  {EIwa;xt,  W:,^)mro)  +  {pf^4>,  <l>t)mn) 

—  {EIWxx,Wxxt)L‘2(rj^)  — (p/ ^t(i,  a:,  0),  Wt)^2(p  )  —  {cBlWxxU'^xxt)L‘^(Va) 

- - - '  ' - _ iJA 

by  integration  by  parts  by  integration  by  parts 


-  {pf'^4>t,  V</>)£,2(n)  +  (/>/  A<j>,  (f>t}L^ci)  -  {pf  X,  0),  wt)i,2(ro) 

'S - ^ - / 

apply  Green’s  formula 

=  {pf'V4>t,  V<j>)xj2(^Q-)  +  {pf  dn<f>,  <f>t)L^(ro)  -  {pf  ^4>t)L^{n) 

-  {pf  X,  0),Wt)L2^r^^  -  {cdIWxxU  Wxxt)L-^(To) 

=  {Pfdn<P,4>t)L^(ro)  -  {pf  a:,  0),  -  {cdIw^xU  '*l^xxt)L'^{To) 

apply  dn4>--wt  on  To 


(11.17) 


-  {pf'(x>t,<i>t)L^ro)  -  {PfM't,X,0),Wt)L2(ro)  -  {cDlWxxt,Wxxt)L^ro) 

—  {^Dl'^xxti'^xx^ L'^{Tff)  • 


Because  —C£)I\\wxxt\Wpt(y^^  <  0,  we  see  that  A  is  in  fact  dissipative. 
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B.  MATHEMATICAL  MODEL  II 


Force 

Figure  3.  2-D  Acoustic  Chamber  for  Model  II 


Our  second  model  is  similar  to  the  first  in  that  we  consider  the  open  loop  case 
in  which  the  wave  equation  is  coupled  with  the  beam  equation  on  one  side  of  a  2-D 
rectangular,  air  filled  cavity  as  shown  in  Figure  3.  Let  fi  =  [0,  a]  x  [0,  /],  5D  =  Fq  U  F 
with  F  =  FiUF2UF3.  Let  the  velocity  potential  in  D  be  given  by  ^  =  <^(t,  x,  y)  and 
the  transverse  beam  displacement  be  denoted  w  =  w{t,x)  as  in  Model  1.  Once  again 
we  take  our  perturbable  boundary  to  be  a  fixed-end  (i.e.,  clamped)  Euler  Bernoulli 
beam  with  Kelvin- Voigt  damping.  However,  instead  of  assuming  hard  wall  boundary 
conditions  along  F  as  we  did  in  the  first  model,  we  now  consider  dissipative  boundary 
conditions  along  F.  Specifically,  the  boundary  conditions  are  taken  to  be 

dn4>  +  a(f>t  =  0  for  {x,  y)  €  F  and  dn(t>  =  Wtit,  x)  for  (x,  0)  €  Fo,  t  >  0 ,  (11.18) 

where  a  is  a  constant  of  proportionality.  As  in  Model  I,  we  take  /(t,  a:)  =  0  for  t  >  0. 
All  other  assumptions  remain  as  stated  in  Model  1.  Hence  the  coupled  system  for 
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Model  II  is  given  by 

(j)tt  =  ,  (x,  y)  e  O ,  f  >  0 , 

dr,<f>  =  -a^t ,  (x,  y)  €  r  ,  t  >  0  , 

dn<t>  =  wt{t,x),  (x,0)  €  To,  t  >  0, 

Pb'u^t-t  +  dl{EIWxx  +  cnlwxxt)  =  -Pf<Pt{t,  x,0) ,  0  <  X  <  a ,  t  >  0,  > 
u;(i, 0)  =  Wa;(i, 0)  =  w{t,a)  =  Wx{t,a)  =  0,  t  >  0, 

<l>{0,x,y)  =  (l>o{x,y)  ,  t0(O,  x)  =  u;o(x) , 

^((0,a:,y)  =  ^^i(x,y)  ,  wt(0,x)  =  wi(x) . 


(11.19) 


The  Hilbert  spaces  H  =  L^{Cl)  x  T^(ro)  and  V  =  x  Ho(ro)  remain 

as  defined  in  the  previous  section. 

In  terms  of  the  sesquilinear  forms  given  by  II.9  and 

Ax(<^,  0  =  /  a  p/  <^e  d'f  for  4>,  ^  €  l'(f2) , 
the  weak  variational  form  of  11.19  is  given  by 

pi{‘<J^tu‘n)  +  +  Piiwir])  =  and  (11.20) 

P2{<f>tt,  0  +  0 P2{4>',  0  —  '^2(wt,0-  (11.21) 


Here  we  pause  to  explain  how  the  boundary  condition  dn<^  +  a4>t  =  0  gives  rise  to  the 
sesquilinear  form  Ai(0,  ^).  Consider  the  wave  equation  =  (PA<f>.  By  multiplying 
through  by  an  appropriately  chosen  trial  function  ^  and  integrating  over  Q,  we  obtain 

I  %4>u^dbj=  /  pjA<i)^doj. 

Ja  Jn 

Looking  just  at  the  right  hand  side  of  this  equation,  we  use  Green’s  formula  to  obtain 


/  Pf^4>^du=f  pf{dncl))(duj  -  [  pjS/(t>-V(duj 
JO,  JdO  Jo 

The  boundary  term  represents 


where  dCl  =  T  U  To  . 


Pf  {dn4>)i  du}  =  pf  wt(  d'f  +  J^pf  (-a)  (j)tC  d'f  =  r2(u;i,  ^)  -  Ai ((?!)(,  ^) . 


17 


The  only  difference  between  the  weak  formulations  for  Models  I  and  II  (equa¬ 
tions  II. 10,  11.11  and  11.20,  11.21  respectively)  is  the  appearance  of  in  equa¬ 

tion  11.21  which,  as  we  have  seen,  arises  because  of  the  damping  along  F.  Just  as  the 
sesquilinear  forms  were  shown  to  satisfy  certain  coercivity  and  continuity  conditions 
in  the  previous  section,  so  too  does  Ai(<^,  ^).  Once  again,  the  Lax-Milgram  Theorem 
assures  us  that  these  weak  formulations  give  rise  to  uniquely  determined  bounded 
linear  operators.  Hence  Model  II  can  be  expressed  as  a  single  first  order  equation 
analogous  to  equation  11.16  in  the  previous  section. 

Letting  $  =  (<^,  w)  and  ^  =  (^,  r;),  such  that  ^  G  V,  and  define  sesquilinear 

forms: 

=  /  p fV (f)  ■  V ^ d(x>  -f  f  Elwxxrjxxd'^  ■, 

Jci  Jtq 

u-2($,'I')  =  /  {cdIwxxIIxx  +  Pf{4>(t,x,0)r]  -  Cit,x,0)w)}  dj  +  /  apf(l>(d'y, 

JTo  Jr 

where  <71,0-2  €  V  xV  — )■  C  (space  of  complex  numbers). 

Once  again,  the  formulations  o-i  and  0-2  satisfy  coercivity  and  continuity  con¬ 
ditions 

5Rcri($,$)  > 

lk.(4’.»)ll  <  o,||4||v||*|k, 

+  Cslklllifr)  =  +  C3||(4|||!(r, , 

Ik2($.«)|l  <  C4||$||v||'»||v. 

Denoting  our  state  variables  by  z{t)  =  (<^,  w)  and  making  use  of  o-i  and  <^2  as 
defined  above,  we  can  express  the  second  order  open  loop  problem  concisely  as 

(Zft(f),  ^)y.,v  4-  cr2{zt{t),  ^)  +  <Ti{z{t),  ^)  =  0  ,  (11.22) 

where  V*  is  the  dual  of  space  V. 

Associated  with  cri  and  <72  are  functionals  Ai ,  A2  such  that  the  weak  and  strong 
formulations  of  the  coupled  system  are  equivalent.  That  is 

(Ai$,  ^)y.,y  =  (7i($,  and  (^2$,  ^)y,y  =  cr2($,  ^) . 
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We  conclude  that  the  system  described  by  equation  11.22  gives  rise  to  an  equivalent 
system  defined  in  terms  of  operators  Ai  and  A2.  That  is 


zu{t)  +  A2Zt{t)  +  Aiz{t)  =  0 .  (11.23) 

The  first  order  system  is  obtained  just  as  it  is  for  Model  I.  We  let  Z{t)  = 
{z(t),zt{t))  and  take  x  G  V  to  obtain 

{Zt{t),x)v\v  =  -CT{Z{t),x)- 

Formally,  this  system  is  equivalent  to  the  strong  formulation 


Z^{t)  =  AZ{t) 


(11.24) 


in  In  matrix  form,  the  strong  formulation  Model  II  is 


<l>t 

Wt 

(t>tt 

.  . 

S.,— V. 


Ut 

where  the  block  matrices  A\,A2  of 


0 

I 

0 

4> 

0 

0 

/ 

W 

0 

0 

0 

4>t 

^dt  -n 

Pb  ^ 

Pb  ^  J 

-  . 

A 

u 

A  are 


(11.25) 


Ai  — 


-c^A  0 

0 

Pb  ^ 


and  A2 


0  0 
n 

Pb  ^ 


Since  energy  is  lost  along  all  four  sides  of  the  cavity  in  Model  II,  it  is  actually 
more  dissipative  than  Model  I,  where  energy  is  lost  only  along  one  side.  Following 
the  procedure  shown  in  11.17  one  finds  that 


{■^Xi  X)u  —  ^Dl\\'>J^xxt\\L‘^(To)  ^{Pj 

where  a  is  the  proportionality  constant  which  appears  in  boundary  condition  11.18. 

Now  that  we  have  variational  formulations  for  Models  I  and  II,  we  turn  our 
attention  to  a  numerical  scheme  often  used  to  obtain  approximate  solutions  to  systems 
of  differential  equations. 
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III.  GALERKIN  APPROXIMATION 

METHOD 


In  order  to  examine  the  exponential  stability  of  approximate  solutions  to  Mod¬ 
els  I  and  II,  we  employ  the  Galerkin  method  to  transform  the  systems  of  partial  dif¬ 
ferential  equations  into  systems  of  ordinary  differential  equations.  A  review  of  this 
popular  computational  technique  follows. 

The  Galerkin  method  is  one  of  the  variational  methods  (e.g.,  Rayleigh- Ritz, 
Least  Squares,  Galerkin,  Collocation,  etc.)  which  all  seek  good  approximate  solutions 
to  many  types  of  boundary-value  and  initial-boundary  value  problems  of  the  form 
Au  =  f  from  a  finite  dimensional  subspace. 

Specifically,  when  given  the  differential  equation  Au  =  /,  where  A  maps  the 
normed  linear  space  X  with  norm  ||  •  Hx  into  the  normed  linear  space  Y  with  norm 
II  •  IIk  and  a  finite  dimensional  subspace  Xn  =  span{4>i,<f>2,  ...,4>n}  of  X,  variational 
methods  seek  functions 


UN  —  Ci<^i  -t-  C2<f>2  +  •••  + 
belonging  to  Xpf  which  minimize 

IIAuiv  -  /||y  +  ||ujv  -  u\\x. 

The  Galerkin  method  is  a  widely  used  technique,  subsuming  both  the  finite  el¬ 
ement  method  and  the  method  of  least  squares.  It  is  worth  noting  that  the  Galerkin 
and  Rayleigh-Ritz  methods  coincide  whenever  the  differential  operator  A  is  linear, 
positive  definite  and  self-adjoint.  However,  there  are  many  important  differential  op¬ 
erators  which  are  either  nonlinear,  not  self-adjoint  or  non-symmetric  for  which  the 
Galerkin  method  is  applicable  while  the  Rayleigh-Ritz  method  fails.  For  example, 
the  Galerkin  scheme  is  applicable  to  many  parabolic  and  hyperbolic  differential  equa¬ 
tions  whereas  the  Rayleigh-Ritz  method  may  not  be  since  the  associated  variational 
problem  may  not  have  a  solution  [Ref.  10]. 
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The  Galerkin  method  yields  a  finite  system  of  equations,  from  a  differential 
equation,  by  discretizing  the  solution  space  rather  than  by  discretizing  the  domain 
and  operators,  as  is  done  in  other  popular  methods.  Approximate  solutions  are  found 
using  a  variant  of  the  method  of  weighted  residuals  (MWR),  where  the  weighting  func¬ 
tions  are  chosen  to  yield  solutions  which  are  a  finite  combination  of  known  functions. 
The  task  at  hand  is  to  determine  approximate  solutions  to  the  linear  (or  nonlinear) 
differential  operator  equation 

Au  =  / 

from  a  finite  dimensional  subspace  of  some  inner  product  space  X  in  which 
the  operator  A  is  defined.  Let  (•,•)  denote  an  inner  product  on  X,  let  Xn  = 
span{^(^,  4>2-,  •••,  =  span{V’('^,  V’iv}  where  X^,  Yn  are  iV-dimen- 

sional  subspaces  of  X.  The  MWR  seeks  an  approximate  solution  ujv  to  Au  =  f  such 
that  UN  G  Xn  satisfies  the  system  of  equations 

{Aun  -  /,<)  =  0,  forj  =  l,2,...,iV.  (III.l) 

Equation  III.l  is  the  inner  product  of  the  residual  {Aun  —  f)  and  the  weighting 
function  integrated  over  the  appropriate  region.  Now  take  un  to  be  a  linear 
combination  of  the  basis  functions  <f>^  such  that 

UN  =  +  C24>2  +  •••  +  (^N<f>N  ■ 

Upon  substituting  this  expression  for  un  into  equation  III.l,  we  find  that  un  must 
satisfy  the  linear  system 

(A(^f ,  t/'f  }ci  =  if,  V’f ) ,  for  i  =  1 , 2, ...,  A" . 

In  the  Galerkin  scheme,  the  weighting  functions  are  taken  to  be  the  basis 
functions  of  the  approximate  solution.  That  is 

V-f  =  for;  =  l,2,....]V. 
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The  approximate  solution  is  obtained  by  solving  the  resulting  system  for  the  coefficient 
functions  Ci{t). 

An  example  will  clarify  the  general  procedure.  Consider  the  following  initial¬ 
boundary  value  problem 


kuxx  =  0 ,  over  F  =  [0, 1] ,  t  >  0 , 

(in.2) 

u(t,  0)  =  0,  1)  =  0, 

(III.3) 

u(0,  x)  =  X  , 

(111.4) 

describing  the  heat  conduction  through  a  thin  one-dimensional  rod  with  no  sources 
of  thermal  energy.  The  temperature  is  held  constant  at  one  end  of  the  rod  (i.e.,  at 
X  =  0);  the  other  end  is  insulated.  Here  we  take  the  thermal  diifusivity  of  the  rod,  k, 
to  equal  one. 

In  operator  notation,  equation  III. 2  can  be  written  as 

Ut  =  Lu 


where  the  operator  L  =  kd^  and,  for  simplicity,  we  take  =  1 . 

Our  first  task  is  to  select  appropriate  basis  functions  which  (i)  possess  the 
smoothness  requirements  of  the  second  order  problem  (i.e.,  <f)^  €  C'^(r))  ,  (ii)  are 
linearly  independent  on  F,  and  (iii)  satisfy  the  boundary  conditions  III.3.  Here  we 
choose 


For  this  example  we  take  N  =  2  obtaining 

=  y-a:  and  <f>^  = 

and  because  we  will  need  the  first  and  second  derivatives  of  these  basis  functions 
later,  we  calculate  them  now: 


=  X -  I ,  =  1 , 

dx4>2  =  -  I,  dl4>2  =  2x . 
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Our  approximate  solution  takes  the  form 


N-2 

u{t,x)  r^UN{t,x)  =  Yj  (^i['t)4>f  {x)  =  Ci{t)<i)^ {x)  C2{t)4>^ {x)  , 

2=1 

where  we  seek  functions  Ci{t).  For  =  2,  we  need  to  select  two  weighting  functions 
ui,u}2  and  introduce  the  spacial  average  (i.e.,  inner  product  or  weighted  integral) 

{u,v)  =  L  uv  d'j , 

such  that 

{u!i,ut)  =  (ui,LuN)  for  z  =  1,2.  (IIL5) 

Note:  The  spacial  superscript  N  is  omitted  below  where  the  meaning  is  clear 
and  an  overdot  denotes  differentiation  in  time. 

For  the  Galerkin  scheme,  the  basis  functions  <l>i  serve  as  the  weighting  func¬ 
tions.  Hence  III. 5  becomes 

un)  =  i4>i,  Lun)  for  f  =  1, 2 


yielding  the  system 


{4>UCi<i>i)  +  {4>l,C2(f>2)  -  -  {cl>i,  C2dl4>2)  =  0, 


{4>2',  Ci(f>i)  -t-  (^2,  ^2^2)  —  Cidl<f>i)  —  {(f>2‘,C2dl<l>2)  —  0, 


or  equivalently, 

{4>2',4>i)  {<t>2i4>2) 


Cl 

0 

{<f>2,dl4>i) 

C2 

0 

(III.6) 


Given  the  particular  boundary  conditions  for  this  problem,  and  because  the  basis 
functions  (jf),  we  selected  satisfy  these  boundary  conditions,  we  find  that  the  inner 
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product  {(t>i,dl(f)j)  is  equivalent  to  the  inner  product  —{dx<i>i,dx4>j)  (i-e.,  L  is  a 
self-adjoint  operator).  To  see  this,  we  integrate  by  parts: 

f  dx  =  (pidxd^jlo  [  dx^i '  dx(j)j  dx  ^ 

Jo  Jo 

where  the  boundary  term  vanishes.  Thus, 

=  -idx<f>i,dx4>j)  ■ 


Although  the  basis  functions  <f)i  were  initially  chosen  so  that  they  belonged  to 
(7^(r),  this  result  indicates  that  less  regular  functions  may  suffice  for  this  problem 
(i.e,  4>i  G  (7^  (r)).  This  fact,  which  is  of  little  use  in  this  simple  example,  is  emphasized 
because  of  its  theoretical  import  as  well  as  its  utility  in  reducing  the  computational 
complexity  of  more  challenging  problems. 

We  now  calculate  the  inner  products  contained  in  equation  III. 6  to  obtain 


2 

15 

6l 

360 

Cl 

+ 

1 

3 

5 

12 

Cl 

0 

61 

68 

C2  _ 

5 

8 

0 

360 

315 

12 

15 

C2 

M  •  A  c  0 


(III.7) 


Thus  the  Galerkin  method  has  reduced  our  problem  of  finding  approximate  solutions 
to  a  second  order  partial  differential  equation  (PDE)  to  that  of  finding  approximate 
solutions  to  a  system  of  first  order  ordinary  differential  equations  (ODEs).  A  much 
simpler  task  indeed!  Multiplying  equation  III. 7  by  M~^  yields 


C\ 

+ 

18.9231 

-5.9077 

Cl 

0 

02 

-12.9231 

7.1077 

C2 

0 

or  equivalently, 


+  18.9231ci  -  5.9077C2  =  0 ,  (III.8) 

b2  -  12.9231ci  d-  7.1077C2  =  0 .  (III.9) 
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To  determine  Ci  and  C2,  we  first  solve  for  in  equation  III. 9  and  then  differentiate 
to  obtain  62  yielding 

C2  =  5^(ci  +  18.9231c,),  (III.IO) 

C2  =  3;^(ci  + 18.92318,).  (III.11) 

Next  substitute  these  equations  for  C2  and  C2  into  equation  III. 9  and  simplify  to  obtain 

C\  T  26.0308ci  T  58.1539ci  =  0  . 

This  is  a  second  order  ODE  with  constant  coefficients  which  is  easily  solved.  We  find 
that  the  general  solutions  for  ci ,  ci  are 

Cl  =  + 

Cl  =  -2.4681  -  23.5628/?e-2^-®®2^% 


where  a,/?  are  constants.  Substituting  these  equations  for  Ci,  Ci  into  equation  III.IO 
we  obtain 

C2  =  2.7853ae-2-^®®“  -  0.7854/3e-^^-^^^^\  (III.12) 

In  order  to  determine  a  and  /?,  we  make  use  of  the  given  initial  condition  given  by 
equation  III.4.  For  (t  =  0,  a;  =  .5)  and  (t  =  0,  x  =  1)  ,  we  know  that  m(0,  .5)  =  .5 
and  u(0, 1)  =  1 ,  respectively.  Hence  we  have 


u(0,.5)  ~  uiv(0,.5)  =  ci(0)^i(.5)  +  C2(0)<?i2(.5)  =  .5 

-1.6515a -0.0151;^  =  .5 

and 

u(0,l)  Ujv(0,l)  =  Ci(0)<^l(l)+C2(0)<^2(l)  =  1 

-2.3569a  +  0.0236/3  =  1. 


From  these  two  equations,  we  determine  that 

a  =  -0.3608  and  (3  =  6.3442 . 


26 


Therefore,  our  approximate  solution  to  equation  III.2  is 

N=2 

UN{t,x)  =  Y,  Ci{t)(t>i{x)  , 

i=l 

where 

Cl  =  +  6.3442e-'*^-^®28i 

C2  =  -1.0049e~2'^®®“-4.9827e"^^-®®2®* 

A  =  T- 

h  =  J-x. 

In  summary,  the  Galerkin  method  simplified  the  task  of  finding  an  approximate 
solution  of  a  PDE  by  recasting  the  problem  as  a  system  of  ODEs  in  a  finite  dimensional 
space.  Because  this  example  was  simply  meant  to  illustrate  the  Galerkin  technique, 
the  solution  obtained  provides  only  a  very  crude  approximate  solution  to  the  example 
problem.  Techniques  for  refinement  of  the  approximate  solution  typically  include 
such  things  as  increasing  the  number  of  basis  functions  and/or  selection  of  different 
basis  functions.  The  interested  reader  is  referred  to  [Ref.  11,  12,  13]  for  additional 
information  regarding  solution  refinement  techniques.  In  Chapter  IV,  we  transform 
Models  I  and  II  into  systems  of  finite  dimension  using  the  Galerkin  method. 
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IV.  FINITE  DIMENSIONAL 
APPROXIMATIONS 


In  Chapter  II,  both  strong  and  weak  variational  forms  were  obtained  for  Mod¬ 
els  I  and  II.  Formally,  at  least,  we  demonstrated  that  these  two  formulations  were 
equivalent.  In  Chapter  III,  we  saw  how  approximate  solutions  to  an  infinite  dimen¬ 
sional  PDE  could  be  obtained  by  employing  the  Galerkin  technique  to  formulate 
the  problem  in  finite  dimensional  spaces.  Key  steps  in  implementing  this  particu¬ 
lar  discretization  technique  included  selection  of  an  appropriately  defined  set  of  basis 
functions,  use  of  these  basis  functions  as  the  weighting  (or  trial)  functions,  calculation 
of  the  inner  product(s)  defined  for  the  discretization  space,  and  finally,  integration 
over  the  spacial  domain  to  transform  the  system  of  PDFs  into  a  system  of  time  de¬ 
pendent  ODEs.  In  this  Chapter  we  extend  these  ideas  to  coupled  systems  of  PDEs, 
to  wit.  Models  I  and  II. 

Our  approach  will  be:  (i)  choose  finite  sets  of  basis  functions  which  span  the 
approximating  solution  spaces,  (ii)  express  the  infinite  dimensional  state  variables 
(u7(t,  x),  <^(t,  a;,  y))  in  terms  of  these  basis  functions,  and  (iii)  use  of  the  weak  varia¬ 
tional  forms  developed  in  Chapter  II  to  obtain  finite  dimensional  representations  of 
Models  I  and  II  necessary  for  our  numerical  work. 

First,  let  denote  the  1-D  basis  functions  which  discretize  the  beam 

and  let  }^j  denote  the  2-D  basis  functions  which  discretize  the  cavity.  For 
the  moment  simply  note  that  there  are  n  -  1  basis  functions  in  and  m  = 

basis  functions  in  where  [rrix  +  l)  and  (rriy  +  l)  represent 

the  number  of  basis  functions  discretizing  the  cavity  along  the  x,  y  axes,  respectively. 
In  Chapter  V,  when  we  look  at  specific  finite  dimensional  approximations  of  Models  I 
and  II,  the  reader  will  better  appreciate  why  these  spaces  have  the  dimensions  given. 

The  basis  sets  {Bf}i~l  and  {5^}^3^span  spaces  and  where  the 
subscripts  b  and  c  denote  ‘beam’  and  ‘cavity’.  That  is  if?  =  span{B?}?-‘  md 
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—  span{B™}^i.  Denoting  the  combined  dimension  of  the  discretized  beam  and 
cavity  spaces  asA^  =  m  +  (n  —  1),  the  approximating  beam  and  cavity  solutions  are 
given  by 

n— 1  m 

w^{t,x)  =  '^w^{t)Bl{x)  and  ^  y). 

7=:1  A:=l 

Notice  that  in  this  form  the  state  variables  are  separated  into  products  of  time  and 
spatial  functions. 

For  application  of  the  Galerkin  scheme,  elements  of  the  ba.sis  sets  and 

serve  as  weighting  functions  for  the  cavity  and  beam,  respectively.  Denoting 
the  product  space  for  the  first  order  system  as  x  ^  restriction  of  the 

infinite  dimensional  first  order  systems  obtained  for  Models  I  and  II  in  Chapter  II  to 
the  space  'H^  x  yields 

<Z'’(().x)«  =  -f(Z''(*),x), 

for  Z^{t)  =  {4’^{t^x,y),w^{t^x)^(l)^{t,x,y),w^{t,x,y)y.  Note  that  (t  is  model  spe¬ 
cific  (Recall  that  we  used  a  to  denote  the  first  order  sesquilinear  vectors  for  both 
Models  I  and  II  in  Chapter  II.).  For  x  =  this  finite  dimensional  first  order 

equation  represents  the  linear  system 

=  (rv.i) 

where 

y^  =  and  (t),- •  ■ 

denotes  the  A/’xl  =  (m-|-n  —  1)  approximate  state  vector.  We  use  an  overdot  to 
denote  differentiation  with  respect  to  time. 

Note:  Below,  and  for  the  remainder  of  this  paper,  the  index  ranges  are  k,i  = 
1 ,  •  •  • ,  m  and  i,j  =  1 ,  •  •  • ,  n  —  1  unless  otherwise  noted. 

Up  to  this  point,  everything  we  have  discussed  in  this  chapter  applies  both 
to  Models  I  and  II.  Now  we  restrict  our  discussion  to  Model  I  for  which  equation 


IV.  1  represents  the  linear  system  shown  below.  The  non-zero  entries  in  and 
represent  block  matrices.  The  rows  of  these  block  matrices  are  generated  by  holding 
£  and  j  fixed  while  varying  k  and  i  as  appropriate  for  the  particular  product  being 
computed.  System  IV.  1  is  given  by 


where  refer  to 

the  sesquilinear  forms  shown  in  II.9.  We  represent  this  linear  system  concisely  as 

0  _  0  Mf 

0  J  [  “  I  -Af  -Af  ^^{t) 


(IV.2) 


with 


0 

0  Mg  _ 


Mg 


'Mg,  o' 

0  Mg,  Y 


II 

1 

o  ' 

II 

1 

o 

- 1 

1 

o 

1 

1 

- 1 
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The  component  matrices  are  given  by 


M^A  =  fvB^-VBTdoj  \  M^A  =[  dlB^dlBJdj 

ik/  L  iij  JTq 

Mii]  =  f  ^B^BTdoj  \  Mii]  =  I  PbBfBJdy 

J  ky£  JCtC  L  liyj  JVq 

1  =  /  Pf^B^-yB^du  \  A^A  =  I  Eld^B^dlBJdj  I  (IV.3) 

Jky£  L  J  2  j  JTo 

'  A^A  =-[  PjBT{^,0)BTdj  \Ag]  =  /  PfBrix,0)BJdj 

Jz/  -'Fo  L  J  kj  -^Fo 

<2  1  =  /  CDldlBfdlBJdj, 

J  z  J  JVq  j 


where  the  finite  dimensional  products  correspond  to  the  infinite  dimensional  sesquilin- 
ear  forms  given  in  II. 9. 

The  finite  dimensional  representation  of  Model  II  is  similar  to  that  given  above 
for  Model  I,  except  the  matrix  contains  the  additional  sub-matrix  which  arises 
because  of  the  boundary  damping  along  F.  For  Model  II,  is  given  by 


^41  ^31 

AN  aN 
^32  ^22  _ 


The  component  matrices  for  Model  II  include  all  of  those  specified  in  IV.3  as  well  as 
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<1  ^UuBfBTd^ 

L  J  k^i 


(IV.4) 


which  corresponds  to  the  weak  form  Ai(^,  ^).  The  linear  system  for  Model  II  is  given 

by 


As  mentioned  previously,  the  rows  of  each  block  sub-matrix  are  generated  by  holding 
i  and  j  fixed  while  varying  k  and  i  as  appropriate  for  the  particular  product  being 
computed. 

The  general  form  of  matrices  and  is  presented  below  to 

help  the  reader  better  conceptualize  the  overall  structure  of  the  system.  The  block 
structure  of  and  is  characteristic  of  the  larger  matrices  and  for 
both  Models  I  and  II.  Also,  the  products  represented  by  (•,  •)  in  the  matrices  below 
correspond  to  those  given  for  Model  I  sub-matrices  A^,  A^,  and  A^  in 

IV. 3.  We  represent  these  matrices  as 
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{BT,Bf) 


{bz,bt) 


{BT,B^_,) 

0 

0 


(  DTi 

0 

0 


{BZ,B^_,) 

0 

0 


(B^B^) 

(B^,B^) 


(-8^1,5?) 


€-iit) 

w^{t) 


(5?,  5^2) 

(5r,5S_i) 


{B-_„B-_,) 


■"’n-2it) 

^n-lb) 


and  = 


(B^,Bf) 

iBf,B^) 


(5^1,52”) 


(B^,B^)  (B^,BT) 


(5^,5? 

(5^,52” 


{B^-i,B^.2) 
(5^-1. 5^-i) 

(5^1.52”) 


<(i) 


(5r,5-_i) 


(5^S^2)  ..-  {B-_„B-_,) 


Dimensions  of  the  matrices  and  sub-matrices  associated  with  Models  I  and  II 


M^,  :  2(m  -|-  (n  —  1))  x  2(m  -(-  (n  —  1)) 

M^,M^,A^,A2  '■  (m -t- (n  —  1))  X  (m -f  (n  —  1)) 


fN  i^N  aN  aN  . 
'115  •‘*^215  ^lU  ^41  • 


m  X  m 


Mg,  Mg,  <,  A"  :  (n  -  1)  X  (n  -  1) 

Ag  :  m  X  (n  —  1)  A^  :  (n  —  1)  x  m 

In  Chapter  V,  we  obtain  specific  numerical  approximations  of  M^  and  for 
Models  I  and  II,  and  examine  the  stability  of  these  finite  systems  to  gain  insight  into 
the  stability  of  the  infinite  dimensional  systems  II.8  and  11.19. 
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V.  SPECIFIC  APPROXIMATIONS  AND 

RESULTS 


In  this  chapter  we  select  specific  basis  functions  for  discretization  of  the  beam 
and  cavity,  discuss  development  of  the  computer  programs  used  to  generate  matrices 
of  the  finite  systems  developed  in  the  previous  chapter,  and  present  results  for  specific 
approximations. 

Cubic  splines  and  tensor  products  of  Legendre  polynomials  are  chosen  as  basis 
functions  for  the  beam  and  cavity  spaces,  respectively  (We  refer  to  the  tensor  prod¬ 
ucts  of  Legendre  polynomials  as  “tensored  Legendre  polynomials”  throughout  this 
paper.).  Since  these  choices  are  by  no  means  the  only  possibilities,  the  interested 
reader  is  referred  to  [Ref.  9]  for  a  discussion  of  alternate  choices  as  well  as  selection 
criteria  which  includes:  smoothness  requirements,  uniform  preservation  of  exponen¬ 
tial  stability  of  approximating  systems,  accuracy,  sparsity  of  system  matrices,  and 
ease  of  implementation. 

The  cubic  splines  used  as  a  basis  for  satisfy  the  smoothness  requirements 
and  are  easily  adapted  to  satisfy  the  clamped  boundary  conditions.  We  construct 
the  set  by  first  partitioning  the  beam  into  n  uniform  intervals  of  step  size 

h  =  Letting  Bf  denote  the  standard  cubic  spline  corresponding  to  this  partition, 
then  the  basis  functions  for  the  beam  discretization  are  taken  to  be 

=  B'^-  2B^  - 

B^  =  Bf  fori  =  2,3,---,n-2 

R:_i  =B:-  2Bl_,  -  2B:^,  , 
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where  the  standard  cubic  splines,  as  defined  in  [Ref.  11],  are  given  by 
{x  -  Ifx  e  [Xi-2,  Xi-i] 

+  ^h'^{x  —  Xi^i)  +  ‘ih{x  —  Xi-iY  —  3(a:  —  if  x  e  [xi_i,Xi] 

Bi{x)  ~  +  Zh?'{xi^i  -  x)  +  3/i(xi+i  -  x)^  —  3(xi+i  —  x)^,  if  x  G  [x^,  x^+i] 

(xi+2  -  xf,  if  X  G  [xi+i,Xi+2] 

0,  otherwise . 

All  of  the  basis  functions  do  in  fact  satisfy  the  clamped  boundary  conditions 

Br(0)  =  d.Bm  =  BrM  =  d^B’(a)  =  0 

for  i  =  1,2,  •  •  •  ,n  —  1  as  can  be  seen  in  Figure  4  below  where  the  interval  [0, 1]  is 
partitioned  into  10  uniform  subintervals.  The  dashed  curves  represent  R”  and 
which  have  compact  support  over  three  intervals.  The  interior  splines,  which  have 
compact  support  over  four  intervals,  are  denoted  by  solid  lines. 

Tensored  Legendre  polynomials  are  used  as  a  basis  for  if™.  As  stated  in  [Ref. 
9],  these  polynomials  produce  smaller,  more  structured  matrices  than  those  obtained 
with  linear  splines  or  finite  elements,  and  the  natural  boundary  conditions  along 
the  cavity  walls  obviate  modification  of  the  basis  elements  to  satisfy  some  essential 
boundary  conditions.  Authors  of  [Ref.  9]  assert,  however,  that  these  benefits  are  not 
as  critical  in  the  2-D  case  as  they  are  in  the  3-D  problem  where  systern  matrices  are 
considerably  larger. 

The  basis  set  of  tensored  Legendre  polynomials  is  obtained  by  forming  the 
product  of  transformed  Legendre  polynomials,  denoted  L^{x)  and  Lj{y),  where  the 
subscripts  i,j  indicate  the  degree  of  the  polynomial,  from  the  interval  [—1,1]  to 
[0,a]  X  [0,/],  respectively.  The  transformed  polynomials  are  obtained  by  substituting 
an  appropriate  linear  transformation  for  x  in  the  recursive  definition  of  the  standard 
Legendre  polynomials: 

Pn+i{x)  =  — +  l)a:P„(x)  -  xPn-i{x)]  (V.l) 

Tt  y "  j. 
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Basis  Splines  (Cubic)  for  the  Beam 


X-axis 


Figure  4.  Cubic  Polynomial  Splines 


with  Po{x)  =  1,  Pi{x)  =  X  . 

For  example,  substituting  ^  for  x  in  definition  V.l  transforms  the  standard 
Legendre  polynomials  from  [-1,1]  to  [0,1].  Orthogonality  of  the  Legendre  polyno¬ 
mials  is  preserved  under  this  linear  transformation  (see  Figure  5).  Recalling  that  the 
set  of  basis  functions  for  the  cavity  is  denoted  we  define  R™  as 

A:  =  1,2,  •  •  •  ,m 

=  for  ^  i  =  0,l,---,m^  +  l  ,  (V.2) 

j  0, 1 ,  •  •  • ,  niy  -j- 1 

where  we  impose  the  condition  i  +j  ^  0  to  eliminate  constant  functions  (i.e.,  exclude 
L^{x)Llj{y)  =  1  for  all  (x,t/))  ensuring  the  set  of  functions  is  suitable  as  a  basis  for 
the  quotient  space.  Hence,  the  dimension  of  the  cavity  space  is  m  =  {rrix  -|-  l)(my  -|- 
1)  —  1.  For  consistency  throughout  this  paper  and  in  our  computational  algorithms, 


Figure  5.  Transformed  Legendre  Polynomials 


the  subscript  k  shown  in  V.2  is  determined  by  holding  j  fixed  while  i  varies.  As  an 
example,  the  indexing  scheme  for  nix  =  2  and  =  2  is 

BT(x,y)  =  L'i(x)VM  Br(x,y)  =  Ll{x)L{(y) 

B2(x,y)  =  Ll(x)L‘„(y)  B^(x,y)- Ll(x)L{(y) 

B^(x,y)  =  L%(x)L[{y)  B?(x,y)  =  L'i{x)L\(y) 

BT(x,y)  =  L\(x)L[{y)  B^(x,y)  =  L%(x)L{(y) . 

Having  selected  basis  functions  for  the  beam  and  cavity  spaces,  we  now  turn 
our  attention  to  the  computation  of  the  various  component  matrices  associated  with 
Models  1  and  11.  All  computations  are  done  using  MATLAB  (MATLAB  is  a  high- 
performance  interactive  software  package  produced  by  The  MathWorks,  Inc.,  for  sci¬ 
entific  and  engineering  numeric  computation.).  Programs  written  for  our  numerical 
work  are  found  in  Appendix  A. 
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Matrices  and  are  all  computed  similarly  since  they  each  involve 

integration  over  the  cavity  space.  In  Chapter  IV  we  defined  these  matrices  as 


^11 


k,i 


J 

=  / 


k,i 


■VBTdu. 


(V.4) 


The  matrices  are  computed  by  a  routine  suggested  in  [Ref.  9].  For  the  moment, 
consider  the  integrand  of  where  R™  =  L“(x)L^(j/)  and  —  Z/“(x)L^(y): 


VBr-VBr  =  (d.BT,d,B'^)-{d.BT,d,BT) 

=  (a.(L‘i').5,(L“L'))  ■  (d,(LlL\),d,{LlL[)) 
=  {L%Lt,  L;d,L‘j)  ■  (L‘,d,L;,  L;a,L',) 

=  (L]L%L‘;d,Ll)  +  {L-Lld,L]d,L[) . 


Because  of  the  structure  of  the  integrand,  we  make  use  of  the  tensor  properties  of  the 
transformed  Legendre  basis  functions  to  construct  M/J,  M^,  and  A^.  Orthogonality 
of  the  transformed  polynomials  reduces  computational  complexity  since 

LjZ/g  =  0  and  J  L“Z/p  =  0  whenever  j^q  and  i  . 

First  construct  fundamental  (m^,  +  1)  x  {rux  +1)  matrices  and  K'^  given  by 

[Mrli,  =  £  and  IKThr  =  £  a.Ll(x)d,Ll(x)  dx 

with  similar  definitions  for  and  (in  fact,  M™  =  Mj”  and  =  R™  when 
rUx  =  ruy  and  [0,  a]  =  [0,/]).  By  orthogonality  of  the  transformed  Legendre  polyno¬ 
mials,  matrices  and  Mp  are  diagonal.  The  matrices  and  are  formed  by 
computing 

=  Mf  (8>  -h  RT  (g)  and  =  . 

The  symbol  (g)  denotes  the  tensor  product,  which  we  accomplish  by  using  MATLAB’s 
kron  function.  The  ordering  in  the  above  definition  is  not  obvious;  however,  attention 
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to  the  indexing  scheme  shown  in  V.3  and  as  described  in  Chapter  IV  delineate  our 
convention.  and  are  obtained  by  removing  the  first  row  and  column  of 
and  reflecting  the  deletion  of  the  constant  function  from  the  basis  set.  Since  we 
take  /?/  to  be  constant  in  this  paper,  is  trivial  to  compute.  Matrices  and 
are  positive  definite  and  symmetric — although  not  sparse.  is  diagonal,  positive 
definite.  The  program  matten.m,  written  to  generate  these  matrices,  computes  the 
transformed  Legendre  polynomials  and  their  derivatives  iteratively.  Integration  of  the 
differentiated  transformed  Legendre  polynomials  is  accomplished  by  using  Gaussian 
quadrature,  while  the  orthogonality  relation  for  Legendre  polynomials,  given  by 


6-a  /-i 


I  ^  Pr{t)P,{t)cH^ 


b  —  a 


2r  +  l 


is  used  to  compute  the  integrals  of  the  translated  Legendre  polynomials  (5rs  denotes 
the  standard  Kronecker  delta:  ^^5  =  0  if  r  ^  5,  =  1  if  r  =  s).  Note  that  is  a 
simple  transformation  factor  which  enables  one  to  use  this  exact  integration  formula 
for  integration  over  an  interval  [a,  6]. 

Since  matrices  M^,  M^,  and  A^,  given  by 

M^A  =  /  dlBfdlBfdj,  \  M^A  =  f  PbBfB'^dj, 

.  Jij  JTo  L  U,J  JTo 

==  Eld^^BfdlBfdj,  CDldlBfdlB]d'^, 

-I  -^ro  L  JiJ  JTo 

all  involve  cubic  spline  functions  or  their  second  derivatives,  they  are  computed  sim¬ 
ilarly.  The  program  myspline.m  is  used  to  generate  the  set  of  basis  splines.  Intrin¬ 
sic  MATLAB  functions  polyder  and  conv  are  used  to  differentiate  and  compute  the 
product  of  the  cubic  splines.  A  simple  program,  polyint.m  ,  performs  the  neces¬ 
sary  integration.  The  programs  used  to  compute  M^,A^,  and  A^  capitalize 
on  the  symmetry  of  the  spline  functions  [matl222.m  computes  A^,  and  A^ ; 
matm22.m  computes  M^).  Because  the  splines  are  equal  to  zero  outside  their  regions 
of  compact  support,  these  four  matrices  become  seven-banded  for  n  >  10.  All  four 
are  symmetric  and  positive  definite  in  construction. 
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The  elements  of  matrices  and  correspond  to  the  integrated  product 
of  the  collapsed  tensored  Legendre  polynomials  (i.e.,  —  L’l{x)V-{y)\y=Q)  and  the 

cubic  splines  over  the  interval  [0,  a].  Recall  these  matrices  are  given  by 

1  =-[  PfBfix)BT{x,0)dj  a,nd  \  Ag]  ^  [  PfBr{x,0)BJ{x)  ,dj. 

Ji/  JTq  L  1  hj  JVq 

These  two  matrices  are  computed  using  the  function  a3132.m  .  Integration,  done  us¬ 
ing  Gaussian  quadrature,  is  accomplished  using  MATLAB’s  intrinsic  function  polyval 
to  evaluate  polynomials  at  translated  Gaussian  knots.  Note  that  computation  of  each 
element  of  Agi  and  actually  requires  three  or  four  integrations  rather  than  just 
one  (three  if  splines  5”“^  or  appear  in  the  integrand;  four  for  integrands  involv¬ 
ing  interior  splines).  This  is  due  to  the  piecewise  construction  of  the  cubic  splines 
over  their  respective  regions  of  compact  support.  For  example,  if  Bf  is  an  interior 
spline,  integration  over  To  is  given  by 

/  p,Bnx)BT(x,<l)dj  =  I  pjBr‘(x)Br(x,0)d7 

Jro  Jto 

PfBf{x)B^{x,0)dx+  r  pfBf{x)B^{x,0)dx 
-2  Jx-1 

fx+1  /•x+2 

A  /  pfBr{x)Bnx,0)dx+  pfBrix)BT{x,0)dx. 

Jx  Jx-\-l 

Although  Agi  and  are  both  full  matrices,  their  computation  is  simplified  since 
each  has  a  well-defined  block  structure.  Further,  since  we  take  pf  to  be  constant,  A^ 
is  precisely  the  negative  transpose  of  A^  (i.e.,  A^  =  (— A^)^). 

The  function  mata41.m  is  used  to  compute 

=  J^apfB^BTdj 

=  apf  Rr(0,  y)BT{0,  y)  dy+  B^{x,  l)BT{x,  1)  dx 
+  1^^  BT{a,y)Br{a,y)dyj  . 

Fundamental  matrices  corresponding  to  integration  across  ri,r2,  and  Fs  are  com¬ 
puted  and  then  summed  to  generate  A^.  Because  of  the  well-defined,  block  diagonal 


41 


structure  of  the  matrices  corresponding  to  integration  across  Fi  and  Fs,  these  ma¬ 
trices  are  computed  simultaneously.  The  matrix  produced  by  integration  across  F2, 
although  not  sparse,  possesses  symmetry  along  diagonals  which  simplifies  its  con¬ 
struction.  Given  the  structures  of  these  fundamental  matrices,  is  symmetric  with 
a  well-defined  block-diagonal  and  symmetric  off-diagonal  construction. 

In  light  of  the  structures  of  the  component  matrices  discussed  above,  we  note 
that  for  both  Models  I  and  II:  (i)  the  matrices  and  are  symmetric  and 
positive  definite  by  construction,  and  (ii)  has  symmetric  and  skew  symmetric 
block  construction  (A  E  is  said  to  be  skew  symmetric  if  A^  =  —A.). 

We  are  now  ready  to  examine  the  stability  of  the  approximation  schemes  devel¬ 
oped  for  Models  I  and  II.  For  our  numerical  work  we  assume  the  following  parameters, 
which  according  to  [Ref.  9],  are  physically  reasonable  for  a  .6m  by  Im  cavity: 

a  =  .6m,  /  =  Im ,  p/ =  1.21^ 

c2  =  117649^,  P6  =  1.355,  EI  =  lZMNm? 

CDi  =  .m^ 

^  see 


Note:  For  Model  II,  we  take  the  proportionality  constant  a  =  1,  where  a  appears  in 
the  boundary  condition  (equation  11.18) 

Qn4>  =  for  (x, y)  €  F ,  t>  0. 


42 


MODEL  I  :  For  n  =  rrix  —  my  =  6,  7,  •  •  • ,  18  the  margins  of  stability  for  the  open 
loop  system  are  listed  in  Table  I.  For  each  n,  the  locations  of  the  eigenvalues,  A, 
are  displayed  in  Figures  6,  7  and  8.  Eigenvalues  were  obtained  using  MATLAB’s 
toolbox  function  eig{A^,M^),  where  matrices  and  are  shown  in  equation 
IV.  1.  Eigenvalues  having  real  parts  with  magnitude  greater  than  1  are  excluded  from 
our  plots  in  order  to  better  see  the  distribution  near  the  imaginary  axis. 

For  small  n,  the  results  obtained  agree  very  favorably  with  those  reported  in 
[Ref.  9].  Comparison  of  the  values  shown  in  Table  I  indicates  that  there  does  not 
appear  to  be  a  uniform  margin  of  stability  between  the  eigenvalues  and  the  imaginary 
axis  (i.e.,  the  data  in  Table  I  does  not  indicate  that  the  maximum  5E(A)  for  the  cases 
tested  is  converging  to  a  limit.).  This  is  what  we  expect  based  on  the  conclusions 
contained  in  [Ref.  5];  however,  we  are  somewhat  hesitant  to  report  this  finding  since 
positive  eigenvalues  appear  in  our  results  for  n  >  17.  These  positive  eigenvalues 
should  not  appear  since  Model  I  is  dissipative,  and  therefore,  all  eigenvalues  of  the 
system  should  lie  in  the  left-half  complex  plane  (i.e.,  9?(A)  <  0).  The  absence  of 
a  clear  margin  of  stability,  as  well  as  the  appearance  of  positive  eigenvalues,  may 
represent  a  numerical/ computational  instability  problem.  We  offer  two  reasons  for 
our  suspicions. 

•  During  the  development  of  the  programs  used  to  compute  the  component 
matrices  M^,A^,  we  were  able  to  delay  the  appearance  of  positive  eigenvalues  by 
incorporating  more  stable  computational  methods.  Our  early  programs  relied  heavily 
on  MATLAB’s  intrinsic  “po/j/-type”  functions  (e.g.,  polyder,  conv,  polyval )  and  our 
simple  polynomial  integration  program  polyint.m.  We  modified  our  code  so  that 
integrations  involving  tensored  Legendre  polynomials — or  derivatives  thereof — is  done 
either  by  using  Gaussian  quadrature  or  by  well-known  properties  of  the  Legendre 
polynomials.  Before  these  changes,  we  observed  positive  eigenvalues  for  n  =  13. 
However,  after  incorporating  these  more  stable  techniques,  positive  eigenvalues  did 
not  appear  until  n  =  17. 
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•  We  first  attempted  to  find  eigenvalues  by  calculating  D  =  and 

by  using  MATLAB’s  eig{D).  However,  for  values  of  n  >  13,  MATLAB  returned  warn¬ 
ings  that  D  was  near  singular.  We  then  sought  to  solve  the  generalized  eigenvalue 
problem  using  MATLAB’s  eig{A^ ,  seeking  a  more  computationally  reliable  al¬ 
gorithm  for  the  problem  at  hand.  Note  that  for  n  <  12,  eig{D)  and  tig{A^ ,M^) 
returned  very  similar  results.  The  consistency  of  the  patterns  shown  in  Figures  6,  7, 
and  8  (all  generated  using  eig  (•,  •)  )  lend  confidence  to  our  belief  that  eig{A^ 
provides  more  reliable  results  than  does  eig  (D)  (Note  that  efy(-,  •)  did  not  return 
any  “near  singular”  warnings  even  when  tested  using  very  poorly  conditioned  ma¬ 
trices.).  Nonetheless,  computational  instability  may  increase  as  n  does  since  MAT¬ 
LAB’s  eig{-,-)  function  uses  a  QZ  algorithm  and,  according  to  [Ref.  14],  some  QZ 
algorithms  destroy  both  the  symmetry  and  positive  definiteness  of  the  semi-definite 
pair  {A^  ^M^). 

Finally,  inspection  of  the  eigenvalue  plots  appearing  in  Figures  6,  7  and  8 
reveals  a  consistency  in  pattern  even  for  n  =  17, 18,  when  positive  eigenvalues  appear. 
Thus  our  computational  approach  does  not  fail  catastrophically  for  a  particular  (large) 
n;  rather,  it  degenerates  as  the  matrix  systems  become  increasingly  singular  as  n 
increases. 
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MODEL  II  :  Results  obtained  for  n  =  =  5, 6,  •  •  • ,  20  are  listed  in  Table  II. 

Eigenvalue  plots  are  displayed  in  Figures  9  and  10.  The  MATLAB  toolbox  function 
eig[A^ ^  M^)  is  again  used.  Eigenvalues  having  real  parts  with  magnitude  greater 
than  30  are  not  displayed  in  order  to  better  see  the  distribution  near  the  imaginary 
axis. 

The  data  reveals  that  while  eigenvalues  lie  further  away  from  the  imaginary 
axis,  as  expected,  given  the  dissipative  boundary  conditions  assumed  along  F  as  well 
as  To,  no  definitive  uniform  margin  of  stability  appears  to  exist.  That  is,  the  values 
shown  in  Table  II  are  creeping  towards  the  imaginary  axis  as  n  increases.  Although 
this  movement  cannot  be  seen  from  Figures  9  and  10,  the  figures  do  reveal  consistency 
in  the  eigenvalue  plots  as  n  increases.  Nonetheless,  this  creeping  phenomena  may  not 
be  an  indication  that  the  infinite  dimensional  system  lacks  uniform  exponential  sta¬ 
bility.  Rather,  the  problem  may  be  related  to  our  numerical/computational  approach 
for  the  reasons  stated  above. 

For  Model  II,  we  see  that  (i)  the  maximum  real  part  of  the  (non-zero)  eigen¬ 
values  lie  further  away  from  the  imaginary  axis  in  this  model  than  they  do  in  Model  I, 
and  (ii)  the  dimension  of  the  approximating  solution  spaces  can  be  increased  without 
the  appearance  of  positive  eigenvalues,  (at  least  \xp  to  n  =  rrix  =  ruy  =  20 — the  extent 
of  our  testing).  Thus  Model  II  is  likely  the  better  choice  for  use  in  formulating  the 
noise  control  problem. 
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n  ^  =  rriy 

max{5R(A)} 

6 

-0.021207 

7 

-0.019998 

8 

-0.015566 

9 

-0.010471 

10 

-0.006223 

11 

-0.004555 

12 

-0.005533 

13 

-0.005440 

14 

-0.006347 

15 

-0.006806 

16 

-0.001130 

17 

+0.050618 

18 

+0.051744 

Table  I.  Model  I:  Margin  between  the  open  loop  eigenvalues  (A)  and  the  imaginary 
axis. 


n  =  rrix  =  rriy 

max{3?(A)} 

5 

-1.218388 

6 

-1.250145 

7 

-1.133252 

8 

-1.388610 

9 

-1.104220 

10 

-1.108291 

11 

-1.088381 

12 

-1.087316 

13 

-1.076611 

14 

-1.072313 

15 

-1.065708 

16 

-1.060883 

17 

-1.055837 

18 

-1.051932 

19 

-1.047767 

20 

-1.044687 

Table  II.  Model  II:  Margin  between  the  open  loop  eigenvalues  (A)  and  the  imaginary 
axis. 
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max{Re(ev)}=-0.006347 


max{Re(ev)}=-0.006806 


Figure  7.  Mod  I:  Eigenvalues  for  n  =  ruy  =  14, 15, 16. 
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Figure  9.  Mod  II:  Eigenvalues  for  n  =  rrix  =  my  =  5, 6, 7, 8, 9, 10, 11, 12. 
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VI 


CONCLUSIONS 


In  this  paper  we  investigate,  by  numerical  approximation,  the  uniform  expo¬ 
nential  stability  of  two  infinite  dimensional  systems  developed  to  model  the  acous¬ 
tic/structure  interaction  of  a  fluid-filled,  rectangular  cavity  (known  to  be  dissipative). 
Model  I  assumes  dissipative  boundary  conditions  along  one  side  of  the  boundary,  while 
Model  II  assumes  dissipation  boundary  conditions  along  all  four  sides  of  the  cavity. 
We  formally  obtain  weak  variational  formulations  for  these  two  models,  express  each 
as  a  finite  dimensional  system  by  discretizing  the  solution  spaces  for  the  acoustic  pres¬ 
sure  (j){t,x^y)  and  transverse  displacement  of  the  beam  w{t,x),  and  use  the  Galerkin 
technique  to  transform  the  systems  of  PDEs  into  systems  of  ODEs.  We  evaluate 
the  uniform  exponential  stability  of  these  systems  by  examining  the  location  of  their 
eigenvalues  in  the  complex  plane.  Eigenvalues  of  these  systems  are  determined  by 
solving  the  generalized  eigenvalue  problem  {\M^  —  A^)y^  =  0.  We  found  that: 


•  The  numerical  approximations  do  not  reflect  the  existence  of  uniform  margins  of 
stability  for  either  model.  The  maximum  real  eigenvalues  do  not  appear  to  be  con¬ 
verging  towards  a  greatest  upper  bound  as  the  dimensions  of  the  finite  systems  in¬ 
crease.  Nonetheless,  our  numerical  results  clearly  indicate  that  Model  II  provides  a 
wider  margin  of  stability  than  does  Model  I  and,  thus,  is  likely  a  better  choice  when 
formulating  the  noise  control  problem. 


•  The  choice  of  cubic  spline  and  tensored  Legendre  polynomials — in  concert  with  the 
use  of  the  Galerkin  method — (i)  simplifies  computation  of  the  component  matrices  of 
and  ^  and  (ii)  contributes  to  the  overall  structure  of  and  .4^  simplifying 
the  computation  of  eigenvalues. 
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Possibilities  for  future  work  to  include: 


•  Investigate  alternate  methods  of  solving  the  generalized  eigenvalue  problem  rather 
than  using  MATLAB’s  eig  Alternate  methods  should  attempt  to  capital¬ 

ize  on  the  structure  of  the  semi-definite  pair  M^).  MATLAB’s  eig  {■,•)  function 
uses  the  QZ  algorithm  and  may  be  destroying  both  the  symmetry  and  positive  defi¬ 
niteness  of  the  pair  {A^ ,  M^)  as  some  QZ  algorithms  do.  An  alternate  approach  to 
the  generalized  eigenvalue  problem  is  suggested  in  [Ref.  14]. 


•  Assume  different  dissipative  boundary  conditions  and  numerically  analyze  the  sta¬ 
bility  of  these  systems  using  various  approximation  schemes.  Consider  models  with 
medium  damping  and/or  different  coupling  mechanisms  between  the  acoustic  and  the 
structure  components  [Ref.  5]. 


•  Investigate  the  numerical  stability  and  preservation  of  exponential  stability  of  the 
approximation  schemes  presented  in  this  paper  with  different  choices  of  basis  functions 
to  discretize  the  beam  and  cavity  solution  spaces,  or  use  different  schemes  altogether. 


•  Investigate  other  mathematical  libraries  such  as  NAG  or  IMSL,  which  may  have 
reliable  subroutines  for  solving  the  generalized  eigenvalue  problem  presented  in  this 
paper. 
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APPENDIX.  MATLAB  FUNCTION  AND 

SCRIPT  FILES 


This  appendix  contains  the  programs  used  to  compute  the  eigenvalues  and 
produce  the  eigenvalue  plots  shown  in  this  paper. 


1.  Below  are  examples  of  the  MATLAB  script  files  which  call  the  various  function 
files  shown  in  this  appendix,  as  well  as  several  intrinsic  MATLAB  files,  for  eigenvalue 
computation  and  plotting. 

Model  I,  Eigenvalue  CompUtatiOB  ***+*=|c***:(::)c  +  ** 
n=6;  mx=n;  my=n;  a=0;  b=.6;  el=l;  rhof=1.21;  rhob=1.35;  c=sqrt (117649) ; 
EI=73.96;  cdl=.001;  [Mll,M21,All]=matten(mx,my,b,el,rhof ,c) ; 

[M12 , A12 , A22] =mat 1222 (a , b , n , El , cdl ) ;  [M22] =matm22 (a , b , n , rhob) ; 

[A31,A32]=a3132(a,b,n,mx,my,rhof) ;  m=(mx+l)*(my+l)-l ;  nml=n-l;  t=m+nml; 
T=zeros(t ,t) ;  M1=[M11  zerosCm.nml) ;  zeros(nml,m)  M12] ; 

M2=[M21  zeros (m,nml) ;  zeros(nml,m)  M22] ;  M=[M1  T;  T  M2]; 

A1=C-A11  zeros (iii,nml)  ;  zeros(nml,m)  -A12]  ; 

A2= [zeros (m,m)  -A31;  -A32  -A22] ;  A=[T  Ml;  A1  A2] ; 

ev6=eig(full(A)  ,full(M))  ;  e6=]iiax(real(ev6))  ;  subplot (2,2, 1)  ;  tt=4*10“4; 

axis([-l  0  -tt  tt]);  hold;  w=[0  0];  q=[-tt  tt]  ;  plot(w,q); 

ww=[-l  0];  qq=[tt  tt] ;  plot(ww,qq);  plot (ev6 , ^ o ' ) ;  title( 'n=mx=my=6’ ) 

4:****************  Model  II,  Eigenvalue  Computation  ♦♦♦♦♦♦♦♦♦sic***** 
n=6;  mx=n;  my=n;  a=0;  b=.6;  el=l;  rhof=1.21;  rhob=1.35;  c=sqrt (117649) ; 
EI=73.96;  cdl=.001;  [Mll,M21,All]=matten(mx,my,b,el,rhof ,c) ; 
[M12,A12,A22]=matl222(a,b,n,EI,cdI) ;  [M22]=matm22(a,b,n,rhob) ; 

[A31  ,A32]=a3132(a,b,n,mx,my ,rhof)  ;  A41]=a41(b,el,iiix,my,rhof  ,c)  ; 
m=(mx+l)*(my+l)-l ;  nml=n-l;  t=m+nml;  T=zeros(t ,t) ; 

M1=[M11  zeros(m,nml) ;  zeros(nml,m)  M12] ; 

M2=[M21  zeros(m,nml) ;  zeros(nml,m)  M22] ;  M=[M1  T;  T  M2]; 

A1=[-A11  zeros(m,nml) ;  zeros(nml,m)  -A12] ;  A2=[A41  -A31;  -A32  -A22] ; 

A=[T  Ml;  A1  A2]  ;  ev6=eig(full(A)  ,full(M))  ;  e6=maLx(real(ev6))  ; 
subplot  (2, 2, 2)  ;  tt=2=t=10''4;  axis([-30  0  -tt  tt]);  hold;  w=[0  0]; 
q=[-tt  tt]  ;  plot(w,q);  ww=[-30  0];  qq=[tt  tt]  ;  plot(ww,qq); 
plot(ev6, 'o') ;  title('n=mx=my=6') 

*  *  *  :|c  4:  *  3(c  5|c  *  %  3|c  3f:  +  +  3f:  3fc  5|e  3|c  +  a|c  +  3jc  +  :|c  %  *  3*:  5jc  *  *  *  3*=  :4: + :ic  3|c  :<c  sj:  4: + +  +  +  :►  +  *  *  :jc  3|c  3*c 
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2.  Function  file  matten.m  computes  and 


***************************************  3|:*************  +  *************3K***3|c* 

function  [Mll,M21,All]=matten(mx,my,b,el,rhof ,c) 

%  [Mil  M21  All]  =  matten(mx,my,b,el,rhof ,c) 

% 

y,  This  function  produces  matrices  Mil,  M21,  and  All — all  are 
%  (mx+l)*(my+l)-l  by  (mx+l)*(my+l)-l 

%  Input :  mx  =  highest  degree  of  Legendre  basis  poly  for  x-axis 
%  my  =  highest  degree  of  Legendre  basis  poly  for  y-axis 
%  b  =  right  end  point  along  x-ax;is  (i.e.,  [0,b]) 

’/o  el  =  right  end  point  along  y-axis  (i.e.,  [0,el]) 
y,  rhof  =  uniform  density  of  fluid 
y,  c  =  speed  of  acoustic  wave  in  fluid 
%  Written  by  Major  J.  M.  Shehaji,  last  update  21  May  95. 

y.  Begin  matten.m 

"It IX  Compute  Mil 
if  mx  ==  my  &  b  ==  el 

mxl=mx+l;  x=ones(l ,mxl) ;  vx=l:2:2*mxl;  intPPx=b*(x./vx) ;  Ma=intPPx; 

%  Determine  Gaussian  weights  (w(i))  k  evaluation  points  (x(i)). 

x=l:l:mx-l;  x=x./sqrt((2*x+l) .*(2*x-l));  j=diag(x,l)+diag(x,-l) ; 

[u  x]=eig(j) ;  x=diag(x);  [x  i]=sort (x) ;  u=u(:,i);  w=u(l,:)."2; 

w=w'.*2;  dx=b/2; 

for  i=l:mx;  s=x(i);  p(i,l)=l,-  p(i,2)=s; 
dpn(i,l)=0,*  dpn(i,2)=l; 
for  j=2:mx 

p(i, j+l)=((2*j-l)*s+p(i, j)-(j-l)*p(i,j-l))/j ; 

dpn(i, j+l)=((2*j-l)*s*dpn(i,j)-(j-l)*dpn(i,j-l)+(2*j-l)*p(i, j))/j ; 
end 

end ;  DPxval=dpn ' ; 
for  i=l:mxl 
for  j=l:mxl 

Ka(i , j)=(2/b)*sum(w^ .*DPxval(i, : ) .*DPxval(j , :)) ; 
end 
end; 

y.y.y.y.y.y.  Mel=Ma  and  Kel=Ka  when  b=el  and  mx=my. 

Mad=diag(Ma) ;  sMad=sparse(Mad) ;  sKa=sparse(Ka)  ; 
MllT=kron(sMad,sKa)+kron(sKa,sMad) ;  [row,col]=size(MllT) ; 


56 


Mll=MllT(2:row,2:col)  ;  All=rhof*sparse(Mll)  ;  */oMll=full(Mll)  ; 

Compute  M21 

M21T=kron(Ma,Ma) ;  t=length(M21T) ; 
M21=sparse(diag((rhof/c~2)*M21T(2:t))) ; 
else 

V/>h  Compute  integrals  of  Legendre  poly’s. 

mxl=mx+l;  x=ones(l,mxl)  ;  vx=l  :2:2=t=mxl ;  intPPx=b* (x . /vx)  ;  Ma=intPPx; 
Mad=diag(Ma) ; 

myl=my+l;  y=ones(l,myl) ;  vy=l :2:2*myl ;  intPPy=el*(y ./vy) ;  Mel=intPPy; 
Meld=diag(Mel)  ; 

y.  Compute  deriv’s  &  eval  ’product’  integrals  of  translated  Legendre’s, 
y,  For  x-axis:  Determine  Gauss  weights  (w(i))  &  eval  points  (x(i)). 
x=l : 1 :mx-l ;  x=x . / sqrt ( (2+x+l) . * (2*x-l) ) ;  j=diag(x , l)+diag(x, -1) ; 

[ux] =eig(j);  x=diag(x) ;  [x  i]=sort (x) ;  u=u(:,i);  w=u(l,:).~2; 

w=w.*2;  dx=b/2; 

for  i=l:mx;  s=x(i);  p(i,l)=l;  p(i,2)=s; 
dpn(i,l)=0;  dpn(i,2)=l; 
for  j=2:mx 

p(i,j+l)  =  ((2*j-l)*s*p(i,j)-(j-l)*p(i,j-l))/j  j 

dpn(i,j+l)  =  ((2*j-l)*s*dpn(i,j)-(j-l)=t:dpn(i,j-l)  +  (2*j-l)*p(i,j))/j  ; 
end 

end ;  DPxval=dpn ’ ; 
for  i=l:mxl 
for  j=l :mxl 

Ka(i,  j)  =  (2/b)*sum(w.*DPxval(i, :)  .*DPxval(j 
end 
end 

y  For  y-axis:  Determine  Gauss  weights  (w(i))  &  eval  points  (y(i)). 
y=l:l:my-l;  y=y ./sqrt ( (2*y+l)  .=t=(2*y-l))  ;  j=diag(y,l)+diag(y,-l)  ; 

[uy] =eig(j);  y=diag(y) ;  [y  i] =sort (y) ;  u=u(:,i); 

wy=u(l,  :)  .■'2;  wy=wy.*2;  dy=el/2; 

for  i=l:my;  sy=y(i) ;  py(i,l)=l;  py(i,2)=sy; 
dpny(i,l)=0;  dpny(i,2)=l; 
for  j=2:my 

py(i, j+l)=((2*j-l)*sy*py(i,j)-(j-l)*py(i,j-l))/j ; 
dpny(i,j+l)=((2*j-l)*sy*dpny(i,j)-(j-l)*dpny(i,j-l)+(2*j-l)*py(i,j))/j; 
end 

end;  DPyval=dpny’ ; 
for  i=l:myl 
for  j=l :myl 

Kel(i, j)=(2/el)*sum(wy.*DPyval(i, :) .*DPyval(j , ; 
end 
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end 

IM  Compute  Mil  &  All 

sMad=sparse(Mad)  ;  sKa=sparse(Ka) ;  sMeld=spajrse(Meld) ; 
sKel=sparsG(Kel) ; 

MllT=kron(sMeld,sKa)+kron(sKel,sMad) ;  [row,col]=size(MllT) ; 
Mll=MllT(2:row,2:col);  All=rhof*Mll ;  Mll=full(Mll)  ; 
y.y.y,  compute  M21 

M21T=kron(Mel,Ma) ;  t=length(M21T) ; 
M21=sparse(diag((rhof/c''2)*M21T(2:t)))  ; 
end  y  End  ’if'  statement, 
end  y  End  matten.m 

*:4c:f:**:|:*3|::|c3|c5f:****+**********************************3|c*3t:**2f:**  +  ****  ********** 
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3.  Function  file  matl222.m  computes  and 


function  [M12,A12,A22]  =  matl222(a,b,n,EI,cdI) 
y,  [M12  A12  A22]  =  ioatl222(a,b,n,EI,cdI) 

y. 

%  Returns  M12,  A12,  and  A22  matrices. 

%  Input:  [a,b]  =  domain; 

y  n  =  no.  of  symmetric  peirtitions  of  interval  [a,b]  ; 

y»  El  =  stiffness  coefficient; 

y.  cdl  =  damping  coefficient. 

y.  Note:  n  >=  4  required. 

y,  Extrinsic  functions  called:  myspline.m 

t 

%  Written  by  Major  J.  M.  Shehaji,  updated  11  April  95. 
y,  Begin  mat  1222. m 

%  Compute  step  size  ’h'  and  generate  'x'  vector. 
h=(b-a)/n;  x=a:h:b; 

%  Compute  the  cubic  spline  basis  set  for  the  beam. 

[B, Bl,Bnml]=my spline (a, b, n) ;  bl=Bl(l,:);  b2=Bl(2,:);  b3=Bl(3,:); 

bnmll=Bnml(l, :) ;  bnml2=Bnml (2, : ) ;  bnml3=Bnml(3, :) ; 
y  Compute  2d  derivative  of  cubic  splines. 

ddbl=polyder(polyder (bl) ) ;  ddb2=polyder (polyder(b2) ) ; 
ddb3=polyder(polyder(b3) ) ;  ddbnmll=polyder(polyder(bnmll)) ; 
ddbnml2=polyder (polyder(bnml2) ) ;  ddbnml3=polyder (polyder (bnml3) ) ; 
[uu  vv]=size(B); 

for  i=13:uu-12  */.  i=13  is  index  of  B(2(l)) 

D2B(i-12, :)=polyder (polyder (B(i, :)))  ; 
end 

yyy  compute  M12  matrix 
if  n  <  4 

errorC'n  >=  4  required') 
el seif  n==4 

mlll=polyint(conv(ddbl,ddbl) ,x(l) ,x(2)) ; 
mll2=polyintCconv(ddb2,ddb2) ,x(2) ,x(3)) ; 

mll3=polyint(conv(ddb3,ddb3) ,x(3) ,x(4)) ;  mll=mlll+mll2+mll3; 

ml21=polyint(conv(ddbl,D2B(l, :)) ,x(l) ,x(2)) ; 
ml22=polyint(conv(ddb2,D2B(2, :)) ,x(2) ,x(3)) ; 
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ml23=polyiiit (conv(ddb3,D2B(3, :))  ,x(3)  ,x(4))  ;  ial2=ml21+ml22+ml23; 

ml31=polyint (conv(ddb2,ddbnmll) ,x(2) ,x(3) ) ; 
ml32=polyint(conv(ddb3,ddbnml2) ,x(3) ,x(4) ) ;  inl3=ml31+ml32; 

m221=polyint(conv(D2B(l, :) ,D2B(l, :)) ,x(l) ,x(2)) ; 
m222=polyint (conv(D2B(2, :) ,D2B(2, :)) ,x(2) ,x(3)) ; 
m223=polyint(conv(D2B(3, :) ,D2B(3, :)) ,x(3) ,x(4)) ; 
m224=polyint (conv(D2B(4, :) ,D2B(4, :)) ,x(4) ,x(5)) ; 
m22=m221+m222+m223+ni224; 

M12=[mll  ml2  ml3;  ml2  m22  ml2;  ml3  ml2  mil]; 
elseif  n==5 

mlll=polyint(conv(ddbl ,ddbl) ,x(l) ,x(2)) ; 
mll2=polyint (conv(ddb2,ddb2) ,x(2) ,x(3) ) ; 

mll3=polyint(conv(ddb3,ddb3) ,x(3) ,x(4)) ;  mll=mlll+mll2+mll3; 

ml21=polyint (conv(ddbl ,D2B(l, : )) ,x(l) 5x(2)) ; 
ml22=polyint (conv(ddb2,D2B(2, : ) ) ,x(2) ,x(3)) ; 

ml23=polyint(conv(ddb3,D2B(3, : )) ,x(3) ,x(4)) ;  ml2=ml21+ml22+ml23; 

ml31=polyint (conv(ddb2,D2B(5, : )) ,x(2) ,x(3)) ; 

ml32=polyint (conv(ddb3,D2B(6, : ) ) ,x(3) ,x(4)) ;  ml3=ml31+ml32; 

ml4=polyint (conv(ddb3jddbnmll) jx(3) ,x(4)) ; 

m221=polyint (conv(D2B(l, :) ,D2B(l, :)) ,x(l) ,x(2)) ; 
m222=polyint (conv(D2B(2, :) ,D2B(2, :)) ,x(2) ,x(3)) ; 
m223=polyint(conv(D2B(3, :) ,D2B(3, :)) ,x(3) ,x(4)) ; 
m224=polyint (conv(D2B(4, ;) jD2B(4, :)) ,x(4) ,x(5)) ; 
m2  2=m2  2 1 +m222+m223+m224 ; 

m231=polyint (conv(D2B(2, :) ,D2B(5, :)) ,x(2) ,x(3)) ; 
m232=polyint(coiiv(D2B(3, :)  ,D2B(6, :))  ,x(3)  ,x(4))  ; 
m233=polyint(conv(D2B(4, :) jD2B(7, :)) ,x(4) ,x(5)) ; 
m23=m231+m232+m233 ; 

m241=polyint(conv(D2B(3, : ) .ddbnmll) ,x(3) ,x(4)) ; 
m242=polyint(conv(D2B(4, :) ,ddbnml2) ,x(4) ,x(5)) ;  m24=m241+m242; 

m331=polyint(conv(D2B(5, :) ,D2B(5, :)) ,x(2) ,x(3)) ; 
m332=polyint(conv(D2B(6, : ) ,D2B(6, : ) ) ,x(3) ,x(4) ) ; 
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m333=polyint(conv(D2B(7, :) ,D2B(7, :)) ,x(4) ,x(5)) ; 
m334=polyint(conv(D2B(8, :) ,D2B(8, ;)) ,x(5) ,x(6)) ; 
m33=m33 1 +m332+m333+m334 ; 

m341=polyint(conv(D2B(6, : ) ,ddbnmll) ,x(3) ,x(4)) ; 
m342=polyint(conv(D2B(7, :) ,ddbnml2) ,x(4) ,x(5)) ; 
m343=polyint(conv(D2B(8, : ) ,ddbnml3) ,x(5) ,x(6)) ; 
ni34=ni341+iu342+in343 ; 

m441=polyint(conv(ddbnmll,ddbnmll) ,x(3) ,x(4)) ; 
m442=polyint(conv(ddbninl2,ddbnml2)  ,x(4)  ,x(5)) ; 
m443=polyint(conv(ddbninl3,ddbnml3)  ,x(5)  ,x(6)); 
m44=m441+m442+in443 ; 

M12=[mll  ml2  ml3  ml4;ml2  m22  m23  m24;  ml3  m23  m33  m34;ml4  m24  m34  m44] ; 
else 

mlll=polyint(conv(ddbl,ddbl) ,x(l) ,x(2)) ; 
mll2=polyint(conv(ddb2,ddb2) ,x(2) ,x(3)) ; 

mll3=polyint(conv(ddb3,ddb3) ,x(3) ,x(4)) ;  mll=mlll+mll2+mll3; 

ml21=polyint(conv(ddbl,D2B(l, :)) ,x(l) ,x(2)) ; 
iiil22=polyint(conv(ddb2,D2B(2, :))  ,x(2)  ,x(3))  ; 

ml23=polyint(coiiv(ddb3,D2B(3, :)) ,x(3) ,x(4)) ;  ml2=ml21+ml22+ml23; 

ml31=polyint(conv(ddb2,D2BC5, :)) ,x(2) ,x(3)) ; 
iiil32=polyint(conv(ddb3,D2B(6, :))  ,x(3)  ,x(4))  ;  ml3=nil31+ml32; 

ml4=polyint (conv(ddb3,D2B(9, :)) ,x(3) ,x(4)) ; 

m221=polyint(conv(D2B(l, :) ,D2B(l, :)) ,x(l) ,x(2)) ; 
m222=polyint(conv(D2B(2, :) ,D2B(2, :)) ,x(2) ,x(3)); 
m223=polyint(conv(D2B(3, : ) ,D2B(3, :)) ,x(3) ,x(4)) ; 
m224=polyint(conv(D2B(4, : ) ,D2B(4, :)) ,x(4) ,x(5)) ; 
m22=m22 1 +m222+m223+m224 ; 

m231=polyint(conv(D2B(2, : ) ,D2B(5, :)) ,x(2) ,x(3)) ; 
m232=polyint(conv(D2B(3, : ) ,D2B(6, :)) ,x(3) ,x(4)) ; 
m233=polyiiit(conv(D2B(4, :)  ,D2B(7,  :))  ,x(4)  ,x(5))  ; 
m23=m231+m232+m233; 

m241=polyint(conv(D2B(3, :) ,D2B(9, :)) ,x(3) ,x(4)) ; 
m242=polyint(conv(D2B(4, : ) ,D2B(10, : ) ) ,x(4) ,x(5)) ;  m24=m241+m242; 
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if  n==6 

m25=0;  %  Since  D2B(13,:)  does  not  exitst  as  defined  above  for  n=6. 

else 

m25=polyint(conv(D2B(4, :) ,D2B(13, :)) ,x(4) ,x(5)) ; 

end 

'/o  Build  M12  matrix: 

wl=[mll  m22*ones(l,n-3)  mil];  w2=[ml2  m23*ones(l,n-4)  ml2] ; 
w3=[ml3  m24*ones(l,n-5)  ml3] ;  w4=[ml4  m25*ones(l ,n-6)  ml4] ; 
M121=sparse(diag(wl)+diag(w2, l)+diag(w3,2)+diag(w4,3)+diag(w2,-l)) ; 
M122=sparse (diag (w3 , -2) +diag (w4 , -3) ) ; 

M12=M121+M122; 

end 

A12=EI*M12;  “/.  Compute  A12 

A22=cdI*M12;  %  Compute  A22 

end  %  End  matl222.m 
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4.  Function  file  matm22.m  computes  . 


function  [M22]  =  matm22(a,b,n,rhob) 

%  [M22]=matm22(a,b,n,rhob) 

•/. 

•/.  The  function  produces  the  (n-l)x(n-l)  M22  matrix  whose  elements 
y.  are  the  integrals  of  rhob*(B(i)*B(j ) )  evaluated  over  the  appro- 
y,  priate  partitions  of  [a,b]  where  i, j=l,2, . .  .  ,n-l.  B  denotes 
%  cubic  splines. 

y. 

y  Input:  a  &  b  =  boundary  of  beam,  [a,b] ; 
y,  n  =  number  of  symmetric  partitions  of  interval  Ca,b]  ; 
y  rhob  =  uniform  mass  density  of  beam 

y  NOTE:  n  must  be  >=  4  for  this  function. 

y. 

y  Extrinsic  functions  called:  myspline.m 

y 

y  Written  by  Major  J.  M.  Shehan,  updated  8  April  95. 
y  Begin  matm22.m 

y  Determine  step  size  and  build  x  vector. 
h=(b-a)/n,*  x=a:h:b; 

y  Compute  cubic  basis  splines  for  beam;  B1=B(1)  &  Bnml=B(n-l) . 

CB,Bl,Bnml]=myspline(a,b,n) ;  bl=Bl(l,:);  b2=Bl(2,:);  b3=Bl(3,:); 
y  Determine  if  'n‘  is  large  enough  and  compute  M22  matrix, 
if  n  <  4 

error ( ’ n  >=  4  required ' ) 
elseif  n==4 

mlll=polyint(conv(bl,bl) ,x(l) ,x(2)) ; 
mll2=polyint(conv(b2,b2) ,x(2) ,x(3)) ; 

mll3=polyint(conv(b3,b3) ,x(3),x(4)) ;  mll=mlll+mll2+mll3; 

ml21=polyint(conv(bl,B(13, :)) ,x(l) ,x(2)) ; 
ml22=polyint(conv(b2,B(14, :)) ,x(2) ,x(3)) ; 

ml23=polyint(conv(b3,B(l5, :)) ,x(3) ,x(4)) ;  ml2=ml21+ml22+ml23; 

ml31=polyint(conv(b2,Bnml(l, :)) ,x(2) ,x(3)) ; 
ml32=polyint(conv(b3,Bnml(2, :)) ,x(3) ,x(4)) ;  ml3=ml31+ml32; 
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m221=polyint(conv(B(13, :) ,B(13, ;)) ,x(l) ,x(2)) ; 
m222=polyint (conv(B(l4, : ) ,B(14, :)) ,x(2) ,x(3)) ; 
m223=polyint(conv(B(l5, :) ,B(15, :)) ,x(3) ,x(4)) ; 
iii224=polyint(conv(B(l6, :)  ,B(16, :))  ,x(4)  ,x(5))  ; 
m22=m221+m222+m223+m224; 

M22=rhob*[mll  ml2  ml3;  ml2  ni22  ml2;  di13  ml2  mil] 
elseif  n==5 

mlll=polyint(conv(bl,bl) ,x(l) ,x(2)) ; 
mll2=polyint (conv(b2,b2) ,x(2) ,x(3)) ; 

mll3=polyint (conv(b3,b3) ,x(3) ,x(4)) ;  mll=mlll+mll2+mll3; 

ml21=polyint(conv(bl,B(13, ;)) ,x(l) ,x(2)) ; 
ml22=polyint(conv(b2,B(14, :)) ,x(2) ,x(3)) ; 

ml23=polyint(conv(b3,B(l5, :)) ,x(3) ,x(4)) ;  ml2=ml21+ml22+ml23; 

ml31=polyint (conv(b2,B(17, :)) ,x(2) ,x(3)) ; 
ml32=polyint(conv(b3,B(18, :)) ,x(3) ,x(4)) ;  ml3=ml31+ml32; 

ml4=polyint(conv(b3,Bnml(l, :)) ,x(3) ,x(4)) ; 

m221=polyint(conv(B(13, :) ,B(13, :)) ,x(l) ,x(2)) ; 
m222=polyint(conv(B(l4, :) ,B(14, ;)) ,x(2) ,x(3)) ; 
m223=polyint(conv(B(15, :) ,B(15, :)) ,x(3) ,x(4)) ; 
m224=polyint(conv(B(l6, :) ,B(16, :)) ,x(4) ,x(5)) ; 
m22=m22  l+m222+m223+ni224 ; 

m231=polyiiit  (conv(B(14, :)  ,B(17, :))  ,x(2)  ,x(3))  ; 
m232=polyint(conv(B(15, :) ,B(18, :)) ,x(3) ,x(4)) j 
m233=polyint(conv(B(16, :) ,B(19, :)) ,x(4) ,x(5)) ; 
m23=m231+m232+m233 ; 

m241=polyint(conv(B(15, :) ,Bnml(l, :)) ,x(3) ,x(4)) ; 
m242=polyint(conv(B(16, :) ,Bnml(2, : ) ) ,x(4) ,x(5) ) ;  m24=m241+m242; 

M22=rhob* [ml 1  ml2  ml3  ml4;ml2  m22  m23  m24;ml3  m23  m22  ml2;ml4  m24  ml2  mil]; 

else 

mlll=polyint (conv(bl ,bl) ,x(l) ,x(2)) ; 
mll2=polyint (conv(b2,b2) ,x(2) ,x(3)) ; 

mll3=polyint(conv(b3,b3) ,x(3) ,x(4)) ;  mll=mlll+mll2+mll3; 


64 


ml21=polyint(conv(bl,B(13, :)) ,x(l) ,x(2)) ; 
ml22=polyint(conv(b2,B(l4, :)) ,x(2) ,xC3)) ; 

ml23=polyint(conv(b3,B(l5, :)) ,x(3) ,x(4)) ;  ml2=ml21+ml22+ml23; 

inl31=polyint(conv(b2,B(17, :))  ,x(2)  ,x(3)) ; 
ml32=polyint(conv(b3,B(18, :))  ,x(3)  ,x(4))  ;  ml3=ml31+nil32; 

inl4=polyint (conv(b3,B(21 , :)) ,x(3) ,x(4)) ; 

m221=polyint(conv(B(13, :) ,B(13, :)) ,x(l) ,x(2)); 
m222=polyint(conv(B(l4, :) ,B(14, :)) ,x(2) ,x(3)) ; 
m223=polyint(conv(B(15, :) ,B(15, :)) ,x(3) ,x(4)) ; 
m224=polyiiit(conv(B(16, :)  ,B(16, :))  ,x(4)  ,x(5))  ; 
m22=in221+m222+m223+m224j 

m231=polyint(conv(B(14, :) ,B(17, :)) ,x(2) ,x(3)) ; 
m232=polyint(conv(B(l5, :) ,B(18, :) ) ,x(3) ,x(4) ) ; 

m233=polyiiit(conv(B(16, :)  ,B(19, :))  ,x(4)  ,x(5))  ;  m23=m231+m232+m233 j 

m241=polyint(conv(B(15, :) ,B(21, :)) ,x(3) ,x(4)) ; 
ia242=polyint(conv(B(l6, :) ,B(22, :)) ,x(4) ,x(5)) ;  m24=m241+m242 ; 

m25=polyint (conv(B(16, :) ,B(25, :)) ,x(4) ,x(5)) ; 

%  Build  M22  matrix: 

wl=[mll  m22*ones(l,n-3)  mil];  w2=[ml2  m23*ones(l,n-4)  ml2] ; 
w3=[ml3  m24*ones(l ,n-5)  ml3] ;  w4=[ml4  m25*ones(l,n-6)  ml4]  ; 
M221=diag(wl)  +  diag(w2,l)  +  diag(w3,2)  +  diag(w4,3); 
M222=diag(w2,-1)  +  diag(w3,-2)  +  diag(w4,-3); 

M22=M221+M222;  M22=rhob*M22; 

end 

end  Vo  End  matm22.m 
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5.  Function  file  aS132.m  conaputes  and 


:|c  :4:  :4:  :|c  :f:  :<c  :f:  **********  +  ************************:*:********  5|c  :f:  3^:  3*:  =|c  **  :f: +  + 3j: 

function  [A31 ,A32]=a3132(a,b,n,mx,my,rhof ) 
y,  [A31  A32]  =  a3132(a,b,n,mx,my,rhof ) 

I 

y  Currently,  this  function  returns  matrices  A31  and  A32.  The  elements 
y  of  these  matrices  correspond  to  the  integrated  product  of  the 
y  collapsed  tensored  Legendre  poly's  (i.e.,  y=0)  and  the  cubic  poly's 
y  over  [a,b] .  Note:  The  cubic  poly's  satisfy  clamped  beam  boundary 
y  conditions.  A32  is  formed  by  computing  -A31'.  Integration  is 
y  accomplished  using  Gaussian  quadrature. 

y 

y  Input:  [a,b]  =  interval  of  integration  (i.e.,  length  of  beam) 
y  n  =  number  of  symmetric  partitions  [a,b]  is  divided  into 

y  mx  =  highest  degree  of  Legendre  poly  in  basis  set  for  beam 

y  my  =  "  "  "  "  "  "  for  cavity 

y  rhof  =  density  of  fluid 

y  Extrensic  functions  called:  legtrans.m  to  compute  Legendre  poly's 
y  myspline.m  to  compute  cubic  splines 

y  Written  by  Major  J.  M.  Shehan,  13  May  95. 

y  Begin  a3132.m 

y  Compute  Gaussian  quadrature  weights  and  knots  for  partitioned  beam. 
h=b/n;  v=0:h:b; 

k=round((4+mx)/2) ;  y  Determine  no.  of  knots. 

x=l:l:k-l;  x=x. /sqrt ( (2*x+l) .*(2*x-l) ) ;  j=diag(x, l)+diag(x,-l) ; 
[ux]=eig(j);  x=diag(x)  ;  [x  i]=sort  (x)  ;  u=u(:,i);  w=u(l , : )  . '"2; 

w=w.*2;  y  'w'  denotes  weights;  'x^  denotes  knots, 
y  Translate  knots  to  appropriate  interval . 
dx=h/2;  u=l:2:2*n;  x=x'; 
for  i=l:n 

X(i, :)=dx*(x+u(i)) ; 

end 

y  Compute  Legendre  &  cubic  poly's. 

[L]=legtrans(b,mx) ;  [B,B1 ,Bnml]=myspline(a,b,n) ; 

y  Evaluate  Legendre's  at  translated  knots  corresp  to  beam  partitions. 

P(1 :n, 1 :k)=ones(n,k) ;  s=0; 
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for  i=2:mx+l 
for  j=l:n 

P(j+n+s, :)=polyval(L(i, :) ,X(j , ; 
end;  s=s+n; 
end 

’/.  Compute  integrals  involving  B(l)  and  B(n-l)  cubic  splines. 
M=zeros(n-l,mx+l) ;  %  Allocate  storage  space. 

ql=w.*polyval(Bl(l, :) ,X(1, :)) ;  q2=w.*polyval(Bl(2, :) ,X(2, :)); 
q3=w.*polyval(Bl(3, :) ,X(3, :)) ; 
rl=w.*polyval(Bnml(l, :) ,X(n-2, :)) ; 
r2=w. *polyval(Bnml(2, :),X(n-l, :)) ; 
r3=w . ♦polyval (Bnml (3 , ; ) , X (n , : ) ) ; 
s=0; 

for  i=l:mx+l 

vl=sum(ql.=i=P(l+s,  :))  ;  v2=sum(q2.*P(2+s, :))  ; 

'  v3=sum(q3.*P(3+s,:)); 

ul=sum(rl.*P(n-2+s, :)) ;  u2=sum(r2.*P(n-l+s, :)) ; 

u3=sum (r3 . ♦P (n+s , : ) ) ; 

M(l,i)=dx*(vl+v2+v3)  ;  M(n-l,i)=dx=t:(ul+u2+u3)  ;  s=s+n; 

end 

’/.  Compute  interior  integrals  (B(i),P(j))  for  i=2,...,n-2  &  for  j=l:mx+l. 
[uu  vv]=size(B); 
s=0;  r=0; 

sl=w.*polyval(B(13, :) ,X(1, :)) ;  s2=w.*polyval(B(l4, :) ,X(2, :)) ; 
s3=w.*polyval(B(15, :) ,X(3, :)) ;  s4=w.*polyval(B(l6, :) ,X(4, :)) ; 
for  i=l : (uu-24)/4 
for  j=l:mx+l 

tl=sum(sl.*P(l+s+r, :)) ;  t2=sum(s2. ♦P(2+s+r, :)) ; 
t3=sum(s3.*P(3+s+r, :)) ;  t4=sum(s4. *P(4+s+r, :)); 

M(i+1 , j )=dx*(tl+t2+t3+t4) ;  s=s+n; 
end;  s=0;  r=l+r; 
end;  '  M=M' ; 

y,  Generate  A31  &  A32  using  M  and  fact  that  P(y=0)  =  (-1)  or  (1)  for  all 
y, collapsed  Legendre  poly's. 

m=(mx+l)*(my+l)-l;  A31=zeros(m+l,n-l) ;  s=0; 

for  i=l:(my+l) 

A31(l+s:mx+l+s,  :)  =  (-l)''(i+l)*M;  s=s+mx+l; 

end 

A31=-rhof*A31(2:m+l, :) ;  A32=-A31’ ; 

end  %  End  a3132.m 
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6.  Function  file  mata41.'m  computes 


3^ :)c  :|c  3jc  3|c  4: 3f:  He  %  3f:  :f:  3)c  4: 9{c  ^  3ic  3)c  :|c  :4c  :i)c  :|c  ^  3(c  :f: :|e  3)c  3|c  3)c  9|c  9ic  3(c ){!:  9jc  3)c  :|c  3{c  ^  3|c  if: 

function  [A41] =mat a41 (b , el , mx , my , rhof ) 
y,  [A41]=mata41(b,el,mx,my,rhof) 

I 

y  Computes  matrix  A41  for  Model  II. 
y  Input  variables:  b  =  right  bdary  of  [0,b] . 
y  el  =  upper  bdary  of  [0,el], 

y  mx  =  highest  degree  of  Legendre  poly  for  beam, 

y  my  =  highest  degree  of  Legendre  poly  for  cavity, 

y  cc  =  proportionality  constant  of  damping  term, 

y  In  this  program,  matrix  A41A  corresponds  to  integration  over 
y  0<=y<=el,  x=0;  A41B  corresponds  to  integration  over  0<=x<=b,  y=el; 
y  and  A41C  corresponds  to  integration  over  0<=y<=el,  x=b. 

y 

y  Written  by  Major  J.  M.  Shehan;  last  update:  23  May  95. 

mxl=mx+l;  myl=my+l;  m=mxl*myl-l;  ml=m+l;  y  Notation  simplification. 

yyy  This  algorithm  cam  be  used  to  computes  A41C. 
y  y=ones(l,myl) ;  v=l :2:2*myl;  intPP=el*(y./v) ;  T=ones(mxl,myl) ; 
y  A41C=2eros(ml,ml) ;  q=0;  g=0; 
y  for  i=l:myl 

y  A41C(l+q:mxl+q, l+g:myl+g)=intPP(i)*T;  q=mxl+q;  g=myl+g; 
y  end;  A41C=2*A41C(2:ml,2:ml) ; 

yyy  This  algorithm  computes  A41A  &  A41C  simultaneously. 
y=ones(l ,myl) ;  vy=l :2 :2*myl ;  intPPy=el*2*(y./vy) ;  A41AC=zeros(ml,ml) ; 
s=0;  q=mxl;  g=myl; 
for  i=l:mxl 
for  j=l :myl 

.Block(i,j)=l+(-l)''(j+l+s)  ; 
end;  s=s+l; 
end 

for  i=2:myl 

A41AC(l+q:mxl+q,l+g:myl+g)=intPPy(i)*Block;  q=mxl+q;  g=myl+g; 
end;  A41AC=A41AC(2:ml,2:ml); 

A41AC(1 :mx, 1 :my)=intPPy(l)*Block(l :mx,2:myl) ; 
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5iXy.  This  alogrithm  computes  A41B. 

x=ones(l,mxl)  ;  vx=l  :2:2*mxl ;  iiitPPx=b* (x . /vx)  ;  tl=intPPx(l) ; 
t2=intPPx(mxl) ;  intPPx(l)=t2;  intPPx(mxl)=tl;  D=diag(intPPx) ; 
A41B=zeros(ml,ml) ;  g=0; 

for  i=l:myl 

A(l:mxl,l+g:mxl+g)=D;  g=g+mxl; 
end 

for  k=l:my 

A41B(k*(mxl)+l: (k+l)*(mxl) , :)=A; 
end;  A41B(l :mxl, : )=A;  A41B=A41B(1 :m, 1 :m) ; 

m  Compute  A41 

A41=-rhof*(spcLrse(A41AC)+sparse(A41B))  ; 
end  y.  End  mata41.m 

3|c  + 3(c  ajc  3|c  *  5^:  3*:  ****  **  :|c  3|c  +  3*6  *  *  3jc  3*c  :|c  *  :|c  *  +*♦*♦*****♦***  +  *******♦*♦******  ♦ 
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7.  Function  file  myspline.m  computes  the  cubic  splines. 


3^:  5*:  :|c  %  5fc  3|c  *  3f:  3jc  :^c  3ic  %  :jc  :ic  5|c  3}c  3|c 3(c  :^c  3}:  5|c  sjc  3ic  3*:  3f:  4: 3^:  3|c  +  3|c  3|c  :f:  ^  +  :jc  :f: 3jc  3jc  :>f:  :|c  3|c * 

function  [B,Bl,Bnnil]  =  myspline(a,b,n) 
y,  [B,Bl,Bnml]  =  myspline(a,b,n) 

I 

'/o  This  function  returns  a  family  of  standard  cubic  splines 
y  defined  over  [a,b] ,  as  well  as  the  "boundary  splines"  which 
y  satisfy  clamped  boundary  conditions.  The  interval  Ca,b]  is 
y  divided  into  'n'  uniform  partitions;  ’h'  is  the  step  size. 

y 

y  Input:  [a,b]  =  interval  along  x-axis; 
y  n  =  no .  of  equispaced  partitions  of  Ca,b] . 

y  Output :  B  =  each  row  of  "matrix"  B  corresponds  to  a  cubic  poly 
y  which  is  defined  only  over  one  step  size  (b-a)/n; 

y  every  4  rows  constitute  a  piecewise  smooth  cubic 

y  polynomial  which  is  non-zero  only  over  4  intervals 

y  (e.g.,  rows  1-4  is  the  first  cubic  poly,  rows  5-8 

y  makes  up  the  second  basis  function,  etc.), 

y  B1  =  left  most  cubic  spline  satisfying  clamped  boundary 

y  conditions. 

y  Bnml  =  right  most  cubic  spline  satisfying  clamped  boundary 

y  conditions. 

y  Written  by  Major  J.  M.  Shehan,  updated  10  April  95. 
y  Begin  myspline.m 
y  Form  standard  cubic  splines. 

h=(b-a)/n;  x=a-3*h:h:b+3*h;  hh=l/h~3;  z=0; 

for  i  =  l:n+3 

B(i+z, :)=hh*[l  -3*x(i)  3*x(i)~2  -x(i)“3]; 

B(i+l+z, : )=hh* [-3  3*h+9*x(i+l)  3*h“2-6*h*x(i+l)-9*x(i+l)"2 

h''3-3*h~2*x(i+l)+3*h*x(i+l)~2+3*x(i+l)''3]  ; 
B(i+2+z, :  )=hh+C3  3*h-9*x(i+3)  9*x(i+3)  “2-6*h*x(i+3)-3*h''2 

h''3+3*h~2*x(i+3)+3*h*x(i+3) '■2-3*x(i+3) “3]  ; 
B(i+3+z,  :)=hh*[-l  3*x(i+4)  -3*x(i+4)‘2  x(i+4)-'3]; 

z=z+3 ; 
end 
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%  Form  exterior  splines  which  satisfy  clcunped  boundary  conditions, 
y.  Form  Bl=B(0)-2*B(l)-2=t=B(-l)  . 

bl=BC7, :)-2*B(10, :)-2*B(4, :) ;  b2=B(8, :)-2*BCll, :) ; 

b3=-2*B(12, :) ; 

Bl=[bl;  b2;  b3] ; 

’/  Form  Bnml  =  B(n)  -  2*B(n+l)  -  2*B(n-l)  . 

[uu  vv]=size(B); 

bnmll=-2*B(uu-ll, :) ;  bnml2=B(uu-7, :)-2*B(uu-10, :) ; 
bnml3=B(uu-6, :)-2*B(uu-9, :)-2*B(uu-3, :) ; 

Bnml=[bnmll;  bnml2;  bnml3] ; 
end  */,  End  myspline.m 

*  :|c  2jc  :<c  :»c  +  *  +  :|c  5tc  3jc  sjc  3(c  3|c  4c  3|c  :4:  +  :fc  +  3|c  *  +  :|c  4:  +  5|c  He  5lc  If:  :|c  *  *  4c  *  5|c  +  4c  sf:  ^  :|c 
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8.  Function  file  legtrans.m  computes  the  transformed  Legendre  polynomials. 


:|c 3^:  if:  3|c  3^:  :|c  3*:  :ic  :4:  +  :f:  3|:  3f:  :Jc  3lc  3jc  ^  :^c  :5|c  :tc  jjc  :f:  :ic  :f:  3|c 3|c  sjc  :jc  5^ 

function  [L]  =  legtrans(b,n) 

‘/o  [L]  =  legtrans(b,n) 

% 

*/,  This  function  produces  a  ^matrix'  L  whose  rows  are  translated 
%  Legendre  polynomials  of  0th  through  nth  degree  defined  on  [0,1]. 
y,  The  0th  degree  polynomial  corresponds  to  the  first  row  of  the 
y  output  matrix,  while  the  nth  degree  polynomial  corresponds  to  the 
y,  n+1  row  (i.e.,  the  last  row)  of  the  output  matrix. 

y. 

*/  Input  arguments:  n  =  highest  degree  of  translated  Legender  poly 
y  desired; 

y  b  =  right  end  point  of  interval  assuming  [0,b] . 

y 

y  Algorithm  written  by  Major  J.  M.  Shehan,  updated  11  April  95. 
y  Begin  legtrans.m 

y  Generate  tranlated  Legendre  polyps  of  0th  &  1st  degree. 

L(l,  :)  =  [zeros(l,n)  1];  L(2,  :)  =  [zeros(l,n-l)  2/b  -1]; 

y  Generate  2d-nth  deg  translated  Legendre  polyps  recursively. 
k=0;  r=n+l; 
for  i=2:n 

d=i+l;  p=[(2*(i-l)+l)*2/b  (2*(i-l)+l)*(-l)] ; 

L(d, :)=(l/i)*([zeros(l,n-i)  conv(p,L(i,n-k:r))]  -  [(i-l)*L(i-l, : )] ) ; 
k=k+l; 
end 

end  y  End  legtrans.m 

3JC  3ic  3j:  *  3|C  :«c  3jc  3|C  3*:  :^c  3j:  3f:  3j:  :jc  3f:  3|C  3^:  sf:  Sf:  *******************  3»:  3f:  3|C  3|c  3f:  3JC  3JC  3f:  3|C  3|C  3f:  3|C  3|C  3*C  3jc  3*:  3|C  3*:  3lc + :j£  3j:  3f: 3fc  :f:  3(C  3jc  3jc  3f: 
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9.  Function  file  polyint.m  performs  polynomial  integration. 


function  [polyint]  =  polyint(p,a,b) 

y.  This  function  integrates  the  polynomial  'p'  over  the  interval  [a,b]  . 
y.  The  polynomial  'p’  is  written  as  a  vector  'v'  with  coefficients  listed 
y,  in  descending  order  (e.g.,  3x“2  +  5x  -  8  ===>  [3  5  8]). 

y. 

y.  Written  by  Major  J.  M.  Shehan  10  Feb  95. 
y,  Begin  polyint.m 

v=[p  0];  y=v./[length(p) :-l:l  1];  polyint  =  polyval(y,b)  -  polyval(y,a) ; 
end  y  End  polyint.m 

*  3jc  :ic  *  *  3)c  3|c  3*c  :|c  *  3lc  *  5|c  3|c  *  4:  *  :f:  3i:  3|c  *  3|c  3|c  ♦  *  *  *  :Jc  *  s(£  :|e  3jc  4:  :|c  *  3|c  3*c  *  3t:  3|c  5|c  5|c  4c  4:  *  4: 4:  :*c  a*:  sf:  3|c  3|c  3|c  3|c  sjc  *  :|c  +  5|c  5jc 
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